NAS A/TM-2015-2 18990 



Entropy Stable Staggered Grid Spectral 
Collocation for the Burgers’ and 
Compressible Navier-Stokes Equations 


Mark H. Carpenter and Matteo Parsani 
Langley Research Center, Hampton, Virginia 


Travis C. Fisher 

Sandia National Laboratories, Albuquerque, New Mexico 


Eric J. Nielsen 

Langley Research Center, Hampton, Virginia 


December 2015 




NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the auspices 
of the Agency Chief Information Officer. It collects, 
organizes, provides for archiving, and disseminates 
NASA’s STI. The NASA STI program provides access 
to the NTRS Registered and its public interface, the 
NASA Technical Reports Server, thus providing one 
of the largest collections of aeronautical and space 
science STI in the world. Results are published in both 
non-NASA channels and by NASA in the NASA STI 
Report Series, which includes the following report 
types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase of 
research that present the results of NASA 
Programs and include extensive data or theoretical 
analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent of 
graphic presentations. 

• TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 
preliminary or of specialized interest, 
e.g., quick release reports, working 

papers, and bibliographies that contain minimal 
annotation. Does not contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 

Collected papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or 

co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. 
English-language translations of foreign 
scientific and technical material pertinent to 
NASA’s mission. 

Specialized services also include organizing 
and publishing research results, distributing 
specialized research announcements and feeds, 
providing information desk and personal search 
support, and enabling data exchange services. 

For more information about the NASA STI program, 
see the following: 

• Access the NASA STI program home page at 
http://www.sti.nasa.gov 

• E-mail your question to help(a>sti. nasa.gov 

• Phone the NASA STI Information Desk at 
757-864-9658 

• Write to: 

NASA STI Information Desk 
Mail Stop 148 

NASA Langley Research Center 
Hampton, VA 23681-2199 



NAS A/TM-2015-2 18990 



Entropy Stable Staggered Grid Spectral 
Collocation for the Burgers’ and 
Compressible Navier-Stokes Equations 


Mark H. Carpenter and Matteo Parsani 
Langley Research Center, Hampton, Virginia 


Travis C. Fisher 

Sandia National Laboratories, Albuquerque, New Mexico 


Eric J. Nielsen 

Langley Research Center, Hampton, Virginia 


National Aeronautics and 
Space Administration 

NASA Langley Research Center 
Hampton, Virginia 23681-2199 


December 2015 




Acknowledgments 


Special thanks are extended to Dr. Mujeeb R. Malik for funding this work as part of 
NASA's Transformational Tools and Technologies (T3) project. This research was also 
supported by an appointment to the NASA Postdoctoral Program at the NASA Langley 
Research Center, administered by Oak Ridge Associated Universities through a contract 
with NASA. Our gratitude also goes to the IT support of our system administrators Steve 
Carrithers and Jim Matthews of the NASA Langley Computational AeroSciences Branch. 


The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 
official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics and 
Space Administration. 


Available from: 


NASA STI Program / Mail Stop 148 
NASA Langley Research Center 
Hampton, VA 23681-2199 
Fax: 757-864-6500 




Abstract 


Staggered grid, entropy stable discontinuous spectral collocation operators of any order are de- 
veloped for Burgers’ and the compressible Navier-Stokes equations on unstructured hexahedral 
elements. This generalization of previous entropy stable spectral collocation work [1,2], extends 
the applicable set of points from tensor product, Legendre-Gauss-Lobatto (LGL) to a combination 
of tensor product Legendre-Gauss (LG) and LGL points. The new semi-discrete operators discretely 
conserve mass, momentum, energy and satisfy a mathematical entropy inequality for both Burg- 
ers’ and the compressible Navier-Stokes equations in three spatial dimensions. They are valid for 
smooth as well as discontinuous flows. The staggered LG and conventional LGL point formulations 
are compared on several challenging test problems. The staggered LG operators are significantly 
more accurate, although more costly to implement. The LG and LGL operators exhibit similar 
robustness, as is demonstrated using test problems known to be problematic for operators that 
lack a nonlinearly stability proof for the compressible Navier-Stokes equations (e.g., discontinuous 
Galerkin, spectral difference, or flux reconstruction operators). 
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1 Introduction 


Beginning with a high-order solution polynomial, numerous approaches exist for constructing dis- 
continuous Galerkin (DG) spectral element methods. Popular variants adopt either the weak (in- 
tegral) or strong (differential) form of the governing equations derived by integrating the equations 
once or twice against a test function. Various interior and interface flux approximations are used 
(e.g., quadrature free fluxes [3], or skew-symmetric operators [4]), as are various quadrature rules 
(e.g., Legendre-Gauss, 1 Legendre-Gauss-Lobatto, or Legendre-Gauss-Radau points). Each design 
choice is motivated by desirable goals, such as efficiency, accuracy, data locality, flexibility, etc. 
Kopriva and Gassner [5] presented a survey of design decisions made when constructing nodal DG 
algorithms, as well as their advantages and disadvantages. 

A popular design philosophy for the incompressible Navier-Stokes equations is to stagger the 
solution and fluxes at independent point sets. For example, Bernardi and Maday use a stag- 
gered approach for the Stokes problem [6]. For the compressible Navier-Stokes equations, Kopriva 
and Kolias [7] collocate the solution variables at the Legendre-Gauss-Chebyshev quadrature points 
(0, • • • , N — 1), and evaluate the fluxes at the Legendre-Gauss-Chebyshev points (0, • • • , N). This 
conservative combination of over-collocated fluxes proves to be more robust in practical problems 
than conventional Legendre-Gauss-Lobatto (LGL) techniques [7]. (The number of flux points ex- 
ceeds the solution points by one, and is reminiscent of the spectral finite volume method of Cai, 
Gottlieb, and Harten [8]). Staggered collocation operators have evolved over the past two decades 
to include a rich set of approaches. For example, the spectral difference (SD) [9, 10] and flux recon- 
struction (FR) [11, 12] approaches, discretize the compressible Navier-Stokes equations in strong 
form on a staggered set of solution (order p ) and flux points (order p+ 1) . The observed convergence 
rate is reported to be (p + 1), for a sequence of nested uniform and quasi-uniform grids. 

An alternate design strategy based on a summation-by-parts (SBP), simultaneous-approximation- 
term (SAT) framework (i.e., SBP-SAT operators), is used in references 1,2 to construct discontin- 
uous collocation spectral element methods of any order, and are referred to as entropy stable (SS) 
discontinuous collocation (DC) algorithms (i.e., SSDC algorithms). Therein, the primary design 
motivation is a semi-discrete spatial operator that supports a nonlinear (entropy) stability proof 
for the compressible three-dimensional (3D) Navier-Stokes equations, on curvilinear unstructured 
hexahedral elements. 

The SSDC operators discretize the governing equations in strong form at the 3D tensor product 
LGL points, and adjoining elements are coupled using a provably nonlinearly stable SAT penalty 
approach technique [13]. The resulting algorithm is similar to the strong form nodal DG method 
reported in reference 14, although it differs in the treatment of the nonlinear Euler fluxes and the 
interface couplings. A novel choice of nonlinear fluxes ensures conservation of mass, momentum and 
energy as well as the entropy inequality within each element; hence element-wise entropy stability. 2 
Carefully constructed interface fluxes then guarantee boundedness of the entropy throughout the 
entire domain. The nonlinear stability is achieved without additional hyper-viscosity dissipation, 
de-aliasing or filtering, and over-integration of the fluxes or solution. Other important differences 
with respect to the DG, SD, and FR schemes appear in the treatment of boundary conditions, which 
are designed to preserve the nonlinear stability of the interior operators (see, references 2, 15). 

1 These points are also referred to as Gauss points in literature. 

2 Dissipation is required for shocked flows to enforce a physical entropy inequality, and density and pressure (or 
temperature) are assumed to remain positive. 
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Entropy stable spectral collocation operators are robust for shocks of moderate strength (e.g., 
normal shock strengths M < 1.75), and are fully consistent with the Lax-Wendroff theorem [16] for 
weak solutions. The robustness stems from a semi-discrete thermodynamic entropy that is provably 
bounded for all time in the L 2 norm, provided that boundary data preserves the entropy estimate of 
the interior operator (and density and temperature remain positive). The nonlinear stability proof 
is sharp; indeed entropy conservative interface fluxes guarantee global entropy conservation (neutral 
nonlinear stability). This sharp estimate is achieved because hyper- viscosity dissipation, de-aliasing 
or filtering of the fluxes/solution is not needed. Furthermore, assumptions of integral exactness 
are unnecessary to justify the proof (commonly used in weak form finite element methods (FEM)), 
because strong conservation form derivatives are approximated using diagonal-norm SBP operators, 
rather than weak form integrals. Thus, over- integration of the nonlinear fluxes is unnecessary to 
more closely approximate integral exactness. 

Although the formulations presented in references 1 , 2 are a significant step towards operational 
entropy stable discontinuous collocation spatial discretizations of any order, noteworthy challenges 
remain: 1) arbitrary collocation points, 2) spatial- ( h ) and order- (p) adaptive refinement of hex- 
ahedral elements, and 3) triangular, prismatic and tetrahedral elements. Herein, an SBP-SAT 
framework is used to develop a staggered grid, entropy stable spectral element formulation that 
includes a broader selection of collocation points. The entropy stable mechanics developed in ref- 
erences 1,2 are extended to include solutions collocated at the Legendre-Gauss (LG) points with 
fluxes computed at the LGL points. The competitive advantages of the new entropy stable stag- 
gered algorithm relative to the algorithms presented in references 1,2 are then established. 

The new staggered operators based on the LG points have several advantages relative to the 
pure LGL operators. First, the integral exactness of the LG points exceeds that of the LGL points: 
(2p + 1) vs. (2 p — 1), and it is shown elsewhere [5] that the LG points have superior accuracy 
properties. Second, the LG points are an interior point distribution and as such, the solution data 
is not duplicated on adjoining element interfaces. Neither are variables collocated at the corners 
of the element. Thus, geometric boundary discontinuities are handled in an integral sense without 
explicit knowledge of the boundary singularity. In more general terms, moving data around an 
element while maintaining entropy stability is an important capability when developing additional 
element types and connectivities. For example, consider the closely related problems of /i-refinement 
at a 2 : 1 element interface compression or a p-refinement interface. These scenarios require data 
movement from adjoining interfaces onto a common intermediate mortar [17]. The quadrature 
points do not in general coincide on either side of the interface. Thus, the nonlinear stability proofs 
presented in references 1,2 do not immediately extend to this extremely important capability. 

Extensive numerical tests presented herein, reveal that the new staggered entropy stable dis- 
continuous collocation (staggered SSDC) operators are significantly more accurate than the LGL 
SSDC operators [1,2], of equivalent polynomial order. They are however, more costly to implement. 
Simple counting arguments based on inviscid and viscous flux evaluations, indicate that the cost of 
the staggered algorithm for a solution polynomial order p is comparable to that of an LGL opera- 
tor [1,2], with a solution polynomial order of (p+1). Preliminary studies using an explicit temporal 
integrator indicate that the increased accuracy of the staggered approach partially offsets the ad- 
ditional cost, particularly at low polynomial order. Further investigation is needed to convincingly 
establish whether over-collocating the fluxes (including the SD [7,9] or FR [11,12] operators) can 
be justified from a cost perspective. An ongoing investigation continues that includes the effects of 
implicit temporal integrators as well as the impact of data movement; computationally intensive 
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yet extremely accurate, low memory footprint algorithms will be competitive in the future. 

The paper is organized as follows. Section 2 includes the general theory of SBP operators, and 
differentiation and interpolation spectral collocation operators. Section 3 includes a brief survey of 
entropy stability at the continuous and semi-discrete level, the data mechanics and the interface 
SAT coupling approach for the staggered SSDC operator in multiple spatial dimensions. Section 4 
includes the mechanics of the staggered SSDC operator. Both energy and entropy proof are pre- 
sented for the ID Burgers’ equation. Section 5 extends the ID staggered SSDC operator to multiple 
spatial dimensions in the context of the compressible Navier-Stokes equations. Section 6 extends 
the staggered SSDC operator to multiple elements. In Section 7, a theoretical cost comparison 
is made between the conventional collocated [1,2] and the staggered SSDC operators. Section 8 
presents numerous numerical studies. The Euler vortex and viscous shock propagation problems 
are used to demonstrate the superior accuracy of the staggered algorithm. The Taylor-Green vortex 
problem and supersonic flow past a 3D square cylinder (Mach = 1.5, and Re = 10 4 ) , are used to 
demonstrate the robustness in the limit of order one discretization errors. Section 9 summarizes the 
results of the paper. Two appendices are included. Appendix A includes a more detailed discus- 
sion on SBP-SAT operators, while a derivation of spectral derivative and interpolation operators 
is included in Appendix B. 

2 Summation-by-parts operators 

2.1 First derivative 

First derivative operators that satisfy the summation- by-parts (SBP) convention, discretely mimic 
the integration-by-parts property 


X L X L 

with cj) an arbitrary scalar test function. At the discrete level, this mimetic property is achieved by 
constructing the first derivative approximation, 770, with an operator in the form 

v = v~ 1 Q, v = v T , C T n>o, c/o, 

Qj = B-Q i B = Diag(— 1,0, ... ,0, 1) , 

where £ is an arbitrary vector. The matrix V can be thought of as a mass matrix (or integrator) 
much like in the finite element framework, or a volume that contains local grid information in the 
context of finite volume or finite difference numerical methods. The nearly skew-symmetric matrix 
Q, is an undivided differencing operator; all rows sum to zero, as do all columns save the first and 
last, which sum to —1 and 1, respectively. The special structure of Q guarantees conservation as 
is proven in the following lemma. 

Lemma 2.1. All differentiation matrices, V, satisfying the SBP convention given in equation (2) 
are discretely conservative in the V-norm. 

Proof. The proof of this lemma can be found in reference 18. □ 
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While the matrix V need not be diagonal, the class of diagonal norm SBP operators play 
a crucial roll in the development of entropy stable (SS) SBP simultaneous-approximation-ternr 
(SAT) operators (see references 1,2,15,19,20). 

The accuracy of the first derivative operator, V, can be expressed as 



!^(x) = D0(x)+T p+ 1 , 

(3) 

where 



x = {xi, 

, ...,X N ), Xl=X L , x N = x H 

(4) 

are the collocated points, and 



0(x) = 

(</>(.ti),</>(x 2 ), • • • ,<Xahv)) T 


d 4> , x 

aL x) = 

( d® , x d<j> # \ T 

(5) 


are the projections of the test function (j) and its derivative onto the grid x. T p+ i is the truncation 
error of the approximation of the first derivative, which is j>th order accurate. Integration in the 
approximation space is conducted using an inner product with the integration weights contained 
in the norm V . 

X H 

J (j)^ dx ~ 4> t VT> q, (6) 

X L 

where 

q(x) = (g(xi),g(x 2 ), ■ • • ,<?(xat)) T , with x = (aq, . . . , x N ) , xi = x L , x N = x H , (7) 

is the projection of continuous variables q onto the grid x. Substituting equation (2) into equation 
(6), the mimetic SBP property is demonstrated, 

= (f) T - Q T ) q = cj) N q N - cf> x q\ - (j) T V T V q. (8) 

2.2 Discontinuous spectral collocation operators 
2.2.1 Differentiation 

Consider the SBP operators constructed at the Legendre-Gauss-Lobatto (LGL) points [21], which 
include the end points of the interval, x L and x H . The complete discretization operator for the 
fourth-order accurate polynomial interpolation ( p = 4) in the standard one-dimensional (ID) ele- 
ment (x L = —1, x H = +1) is illustrated in Figure 1. In this figure, the solution points are identified 
with • and the flux points are identified with j . The latter points are similar in nature to the control 
volume edges employed in the finite volume method and are used to prove the nonlinear stability 
(entropy stability) as briefly shown in Section 3.3, (see references 1,2 for a more detailed discussion). 


Define the Lagrange basis polynomials relative to the N discrete LGL points, x, as 

rN X — Xk 


%w=n la 


k^j x j x k 


1 < j < N. 


(9) 


7 


h 

U 3 


h 


fo h h 

h h 


fi /s 
h /s 


u 1 


U 2 


U4 


U 5 








®1 


X2 


£3 


X4 


^5 

S 0 

Xi 
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X 3 


X4 

x 5 
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10 
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45 
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Figure 1. The one-dimensional discretization for the fourth-order accurate polynomial interpolation 
( p = 4) LGL collocation is illustrated. Solution points are identified with • and the flux points are 
identified with I. 


Assume that a smooth and (infinitely) differentiable function /(x) is defined on the interval 
x L = — 1 < x < 1 = x H . Reading the function / and its derivative ^ at the discrete points, x, 
yields the vectors 


f ( x ) = (/(^t), f{x 2 ), • • • , f(x N -i),f(x N )) 


df 


dx 


df 


df 


df 


df 


t:( x ) = iz( x 1 )’ ^(* 2 ), ■ ■ ■ , 1 ~{xn- 1 ), -t-(x n ) 


dx 


dx 


dx 


dx 


T 


(10) 


The interpolation polynomial /at(x) (of order p = N — 1) that collocates f{x) at the discrete 
points, x, is given by the contraction 


f{x) « f N (x) = L(x; x) T f(x), 


( 11 ) 


where L(x; x) is a column vector whose components are the Lagrange basis polynomials relative to 
the nodes x (i.e. , Lj(x) in equation (9)). Note that the explicit dependence of L on the independent 
variable x and the set of point x is indicated for completeness. 3 

Theorem 2.1. The derivative operator that exactly differentiates an arbitrary p-th order polynomial 
(p = N — 1 ) at the collocation points, x, is 

v = (iij) = , ( 12 ) 

where ^f(xj) denotes the derivatives of the Lj Lagrange basis polynomial with respect to x evaluated 
at the collocated node Xj. This element corresponds to the element in the j-th column and i-th row 
of the differentiation matrix V. 

Proof. The proof to this theorem can be found in Appendix B (see also reference 21). □ 


A representation of the differentiation operator T>, which satisfies all the requirements for being 
an SBP operator is given in the following theorem 


3 /jv is a polynomial of order p in the independent variable x. 

8 


Theorem 2.2. The derivative operator that exactly differentiates an arbitrary p-th order polynomial 
(p = N — 1) at the collocation points, x, can be expressed as 

V = V^ 1 Q (13) 

with 

'P = E£ l (^; x ) l (^; x ) t ^ , Q = E^ l ( 7 7«; x )^( j ?/; x ) t ^) ( 14 ) 

where rp and cue, 1 < l < N, are the abscissae of the LGL points and their quadrature weights, 
respectively. L(x; x) is a column vector whose components are the Lagrange basis polynomials 
relative to the discrete nodes x (i.e., Lj(x) in equation (9)) 

Proof. The proof to this theorem can be found in Appendix B. □ 


The matrix V in (14) is symmetric and positive definite for any vector x [21]. When the 
LGL points are used for x, then V is a diagonal approximation (i.e., the so-called “mass lumped” 
approximation) of the full 'P-norrn (see Appendix B). A diagonal norm SBP operator is necessary 
to achieve strict entropy stability [22,23]. This constraint on the norm V is reiterated in section 
3.3.1. 


2.2.2 Interpolation 

Define on the interval — 1 < x < 1, the vectors of discrete point, 

X = (21,22, • • • ,2M-1,2 M ) T , “I < Xi,X 2 , ■ ■ 
X = (21, 2 2 , • • • , 2at_i, 2at) T , -1 < Xi, 2 2 , ■ ■ 


,xm-i,xm < l ; 

■ , XN-l, Xn < 1. 


(15) 


Herein, the discrete points x and x are the Legendre-Gauss (LG) points (i.e., the so called Gauss 
points) and the LGL points, respectively. All the scalars, vectors, and matrices associated to the 
LG points are denoted with a “tilde” symbol. Next, define the interpolation operators that move 
data between x and x: 

ZlGL^LG = V 1 P- LG- LGL, 

>T 

■'LG— LGL, 


’Plg^lgl = V L Pi 


(16) 


V Plgl- 


• LG 


= 1, 


■T 

LG^LGL 


V, 


where 

l 

Plg-lgl = J L( 2 ; x) L( 2 ; x) t dx. (17) 

-l 

In Appendix B.3, we prove that these polynomial interpolation operators exist and satisfy the 
relations (16), provided that the LGL points are of higher polynomial orders than the LG points 
(i.e., N > M). 
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3 Entropy consistent and entropy stable SBP operators 

3.1 Governing equations 

Consider the three-dimensional (3D) compressible Navier-Stokes equations for a calorically perfect 
gas expressed in the form 


dq dfi _ df[ V) 

dt dxi dxi 


iGO, 


Bq = g^ B \x,t), x e <90, 
q(x,0) = g (0 \x), ieSl, 


t e [o, oo), 
t e [o, oo), 


(18) 


where the Cartesian coordinates, x = (aq, X 2 , xfi) , and time, t, are independent variables, and 
index sums are implied. The vectors q, fi , and /• ' ^ are the conserved variables, and the conserved 
inviscid and viscous fluxes, respectively. Without loss of generality, a 3D box 

n = [xf,xf] x [x%,x$] x 


is chosen as our computational domain with <90 representing the boundary of the domain. The 
boundary vector g ^ is assumed to contain linearly well-posed Dirichlet and/or Neumann data. 
Herein, we have omitted a detailed description of the 3D compressible Navier-Stokes equations 
because it can easily be found in literature. 


3.2 Continuous analysis 

Consider the (nonlinear) compressible Navier-Stokes equations given in equation (18). This system 
of incomplete parabolic partial differential equations (PDEs) have a quadratic or otherwise convex 
extension of its original form, that when integrated over the physical domain, P, depends only 
on boundary data and dissipative terms. This convex extension yields the entropy function and 
provides a mechanism for proving the stability in the L 2 norm of the nonlinear system of PDEs 
(18). In fact, Dafermos [24] showed that if a system of conservation laws is endowed with a convex 
entropy function, S = S(q), a bound on the global estimate of S can be converted into an a priori 
estimate on the solution vector q (e.g., the solution of system (18)). 

Definition 3.1. A scalar function S = S(q) is an entropy function of system (18) if it satisfies 
the following conditions: 

• Differentiation of the convex function S(q), simultaneously contracts all the inviscid spatial 
fluxes as follows 


di dfi = dsdfi dq^ = m dq_ = m . = . 

dq dxi dq dq dxi dq dxi dxi ’ % ’ ’ 

The components of the contracting vector, dS/dq, are the entropy variables denoted as w T = 
dS/dq. Fflq) are the entropy fluxes in the i-direction. 
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• The entropy variables, w, symmetrize system (18) if w assumes the role of a new independent 
variable (i.e., q = q(w)). Expressing equations (18) in terms of w yields 

dq dfi df j V ^ dq dw df t dw d ( ^ dw \ . ^ 

dt + d^i ~ dxt = dw ~dt + Ihv dxl ~ d^j)~ ’ ? = ’ ’ ’ 

with the symmetry conditions: dq/dw = (dq/drv) T , dfi/dw = ( dfi/dw) T andcij = c-. 

Because the entropy is convex, the Hessian d 2 S/dq 2 = dw/dq is positive definite 4 , 

ff 2 9 

C T ^C>o, vc/o, ( 21 ) 

and yields a one-to-one mapping from conservation variables, q, to entropy variables, w. Likewise, 
dw/dq is symmetric positive definite (SPD) because dq/dw = (dw/dq)- 1 and SPD matrices are 
invertible. The entropy and corresponding entropy flux are often denoted an entropy- entropy flux 
pair, (S, F) [25]. 

The symmetry of the matrices dq/dw and dfi/dw, indicates that the conservation variables, q, 
and inviscid fluxes, ft, are Jacobians of scalar functions with respect to the entropy variables, 

T _ dip T _ dfi 
y. 5 J i o’ 
ow aw 

where the nonlinear function, ip, is called the potential and ifi are called the potential fluxes [25]. The 
potential and the corresponding potential flux are denoted a potential-potential flux pair, (ip,f>) [25]. 

Just as the entropy function is convex with respect to the conservative variables ( d 2 S/dq 2 is 
SPD), the potential function is convex with respect to the entropy variables, w. 




The two elements of Definition 3.1 are closely related, as is shown by Godunov [26] and Mock 
[27]. Godunov proves that: 


Theorem 3.1. If equation (18) can be symmetrized by introducing new variables w, and q is a 
convex function of <p, then an entropy function S = S(q) is given by 

p = w T q — S, (23) 


and the entropy fluxes Ffiq ) satisfy 


A = w T fi - hi. 


Proof. The proof of this theorem can be found in references 26,28. 


(24) 

□ 


Mock proves the converse to be true: 

Theorem 3.2. If S = S(q) is an entropy function of system (18), then w T = dS/dq symmetrizes 
(18). 

Proof. The proof of this theorem can be found in references 27,28. □ 

4 Thc Hessian, d 2 S/dq 2 , is actually symmetric positive definite (SPD). 
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Contracting system (18) with the entropy variables, w, results in the differential form of the 
(scalar) entropy equation, 


dSdq dSdfr = ds m = d^df^ = JL ( t av)\ _ (dw\ T av) 

dq dt dq dx; dt dxi dq dxi dxi V * / \ dx t ) 1 



(25) 


Integrating equation (25) over the domain yields a global conservation statement for the entropy 
in the domain 


d_ 

dt 



T AV) 
w fl 



v>l«j w Xj dx. 


(26) 


n 


n 


References 22,23 prove that the five-by-five matrices c t] in the last term in the integral are positive 
semi-definite. Note that the entropy can only increase in the domain based on data that convects 
and diffuses through the boundaries, dll. The sign of the entropy change from viscous dissipation 
is always negative. 


Remark 3.1. The scalar equation (26) is an integral statement of entropy conservation hut it is 
not strictly valid in the presence of discontinuities (i.e., shocks). In fact, it does not account for 
the dissipation of the entropy at a shock. Although the precise amount of entropy dissipated at a 
discontinuity is not known a priori, what is known is the sign of the jump in entropy. Thus, a 
general but not sharp statement of the global behavior of entropy in the entire domain is 


d_ 

dt 


J 


5dx < 


T AV) 


W l 



n 


(27) 


Remark 3.2. A sufficient condition to ensure the convexity of the entropy function S = S(q) (and, 
hence, a one-to-one mapping between the entropy variables, w, and the conservative variables, q, is 
that p, T > 0 (for the proof see for instance, Appendix B.l in reference 19). Expressly: 

C T ffc T >0, VC f o, p,T> 0. 

dq- 

This (physical and mathematical) restriction on density, p, and temperature, T, weakens the en- 
tropy proof, making it less than full measure of nonlinear stability. Another mechanism must be 
employed to bound p and T away from zero to guarantee positivity; positivity preservation will not 
be considered herein. 


3.3 Overview of the semi-discrete entropy analysis for the LGL points 

The semi-discrete entropy estimate is achieved by mimicking term by term the continuous estimate 
given in equation (26). As for the continuous case, the nonlinear stability (entropy stability) 
analysis begins by contracting the discrete entropy variables, w T , with the semi-discrete version of 
the system (18) (see for instance, references 1,2). (For clarity of presentation, but without loss of 
generality, the derivation is simplified to one spatial dimension. Tensor product algebra allows the 
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results to be extended directly to three-dimensions.) The resulting global equation that governs 
the senri-discrete decay of entropy is given by [1] 

w T Vq t + wAf = w T Af (V) + w T g (B) + w T g (/nt) , (28) 

where 

w = (w(q 1 ) T ,w{q 2 ) T ,...,w(q N ) T ^j . 

The source terms g ^ and gd nt ) contain the enforcement of boundary and interface conditions, 
respectively. (Herein, the solution between adjoining elements or cells is allowed to be discontinuous. 
Therefore, interface penalties g 1 Int ' 1 are needed to patch interfaces together). The entropy variables, 
w, are defined at the solution points whereas the quantities with an over-bar, i.e., f and F, are 
defined at the flux points (see Figure 1). The analysis of each senri-discrete term is presented 
elsewhere [19]. The stability of the viscous ^w T Af^' ^ terms follow immediately by conventional 
approaches for diagonal norm SBP operators and will not be repeated herein. The analysis of 
the time derivative and inviscid terms is summarized next, while the form of the penalty terms is 
presented in Section 6. 


3.3.1 Time derivative 


The time derivative in (28) is in mimetic form for diagonal norm SBP operators. The entropy 
variables are defined by the expression w T = dS/dq. Define the diagonal matrices <9S/<9q = W = 
Diag[w]. Since V is a diagonal matrix (see equation (14)) and arbitrary diagonal matrices commute, 
the senri-discrete rate of change of entropy becomes 


= UwP^ = = Up**! = Up* 


dt 


dt 


dt 


dq dt 


dt 


where 

1 = ( 1 , 1 ,..., 1 ) T , 

is a vector with N elements. 5 


3.3.2 Entropy consistent inviscid fluxes 

The inviscid portion of equation (28) is entropy conservative if it satisfies 

w T Af = F(q N ) - F( qi ) = F(q N ) - F( qi ) = 1 T AF. (29) 

Tadnror [25] uses the following rational to construct entropy conservative (entropy consistent) 
three-point centered operators satisfying equation (29). Consider the term Wi(f i — / i _ 1 ) in the 
vector relation (29), which can be rewritten as [w T Af— AF] = 0. Adding and subtracting equivalent 
terms yields the expression 

Wi (7 i ~ 7i-i) = [| («>i+ 1 + w%) l_i ~ \ (m + Wi- 1) /i-i] /o n x 

- [\ (wi + 1 - u>i) fi + 3 (wi - Wi- 1 ) fi_i) . 

5 N is the size of the ID grid x. 
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Note that in this section the subscripts i — 1, i and i + 1 are used to denote a scalar or vector 
quantity at the i — 1, i or i + 1 collocated point, and should not be confused with the subscript 
used, for instance, in (18) to indicate the coordinate direction. The first two bracketed terms on 
the right-hand-side of equation (30) telescope across the domain in their current form (i.e., their 
contributions within an element sum to zero). The last two terms telescope provided they satisfy 
a consistency condition of the form 

(m + 1 - w % ) 7i S) = A~ ipi-i, 2 < i < N (31) 

—(S) 

(modulo slight changes at the end points of the domain) , where f i denotes an entropy conservative 
(or entropy consistent) flux. 



A general strategy for constructing an entropy conservative flux, f i , that satisfies the point- 
wise conditions 

(w i+ 1 - Wi) f\ S) = A + 1 ~A, i = 1, 2, . . . , N - 1 ; A = A, ipN = ipN (32) 

—(S) 

is presented elsewhere [23]. Herein, the flux f i is based on linear combinations of qij- weighted, 
two-point entropy conservative fluxes f s = f s (ii£,Uk), which satisfy the following relation: 

(we ~ w k ) Js ( u e , uk) =A- V’fc- (33) 


The dyadic shuffle conditions given by equation (33) are known to exist for Burgers’ equation and 
the Euler equations [25,29]. 

The following theorem summarizes the work given in reference 23, and provides the general 
—(S) 

formula for constructing f i of any order from a linear combination of dyadic entropy conservative 
fluxes / s(ue,u k ). 

Theorem 3.3. A two-point high-order accurate entropy conservative flux satisfying equation (32) 
with formal boundary closures can be constructed as 

N i 

7i S) = £ 2 7s («*«*), l<i< AT-1, 

k=i + 1 e=\ 


where f s (ug,Uk) is any two-point non- dissipative flux function that satisfies the entropy conserva- 
tion condition given by equation (33). The two-point high-order accurate entropy conservative flux, 

_(g) 

f i , satisfies an additional local entropy conservation property, 


w T 'p- 1 Af (S) = V- 1 AF = ^ (q) + T p+ 1 , 


or equivalently, 
where 


Wi 


• T (fi S) - /S) = (Fi - Fi-ij , 1 <i<N, 


N 


Fi = EE qtk \ (wg + w k ) T f s (u e , u k ) - (A + A) 

k=i+ 1 1= 1 


1 < i < N - 1. 

Proof. For brevity, the proof is not included herein, but is reported elsewhere [23]. 


(34) 

(35) 

(36) 

□ 
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4 Stability on staggered grids: Burgers’ equation 


An energy /entropy analysis of the ID Burgers’ equation is presented before that of the compressible 
Navier-Stokes equations. 

4.1 Data mechanics of the staggered grid approach 

Define a staggered grid algorithm (Sta-Grd-Alg) for building discrete differentiation operators using 
two sets of collocation points: x and x of dimension M and N, respectively (see Figure 2). The 


proposed algorithm is 

similar to that proposed by Kopriva and Kolias [7]. 

Assume that the time- 

dependent solution is 

stored at the points 

x. Furthermore, 

assume 

that the extrema of x coincide 

with the endpoints of the domain: 

Xl = 

x L , x^v = x H , 

to facilitate imposition of interface or 

boundary data. 







/o fl 


/a 


fa 
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Figure 2. The one-dimensional discretization for the fourth-order accurate polynomial interpolation 
(p = 4) with the staggered approach is illustrated. Solution points x are identified with x and 
auxiliary points x are identified with •. Flux points x (used to prove the entropy stability) are 
identified with |. 

Discrete differentiation of first or second order spatial terms (e.g., df /dx or df^/dx), by using 
the Sta-Grd-Alg is accomplished as follows: 

• Interpolate the discrete entropy variables from x to x. 

• Build the nonlinear fluxes / and ) on the set of points x. 

• Build the interface and/or boundary penalties at the extrema of x. 

• Differentiate the fluxes on x, and impose the penalties by using the SAT approach. 

• Interpolate the discrete flux derivatives and penalties back to x. 

• Advance the solution with a time integration scheme by using the interpolated flux derivative 
on x. 

Tensor product arithmetic extends the approach directly to three spatial dimensions. 6 An SBP- 
SAT stability proof is now presented for the Sta-Grd-Alg. It is valid for all diagonal-norm SBP 
operators. 

G The Sta-Grd-Alg is valid for other grid distributions that do not support tensor product arithmetic. 
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Define “tilde” variables and operators that act on the set of points x (e.g., u, V and V ), and the 
pair of interpolation operators Ilg^lgl and Ilgl^lg introduced in section 2.2.2. The Ilg^lgl 
operator transfers data from x to x, while the Ilgl^lg transfers data from x to x. Define the 
interpolated solution vector u, and the diagonal velocity and viscosity matrices [n] and [e] as 

u = Ilg^lgl u ; [n] = Diag[2L G ^LGLu] ; [e] = Di&g[l LG ^LGiA\ ; (37) 

and the boundary operator nomenclature 

e (x L ) = (1,0, ••• ,0)1, 
u (x L ) = uL=_i = u T e (x L ) , 

V ' ' V -r r (38) 

{Vu) ( x L ) = (Du) | x= _i = (Du) T e ( x L ) , 

g ( x L ) = g|*=_i = g T e ( x L ) , 

with similar definitions for e (x H ) , u ( x H ) , (Du) {x H ) , and g ( x H ) . The subscript M in (38) denotes 
the size of the vector e (x L ) . 


4.2 Stability of Burgers’ equation 


Conventional energy estimates as well as entropy analysis exists for Burgers’ equation for all diag- 
onal norm SBP operators [14, 23, 25] . Comparison of the two approaches provides insight on how 
to proceed with the analysis of the compressible Navier-Stokes equations. 

Consider the ID viscous Burgers’ equation 


. BJin) = 0(ef iv) (m ,,, = t 

dt dx dx 2 



du 
dx ’ 


x e [ x L , x H ], 


t € [0, oo), 


u ( x L ,t ) + | u (x L ,t)\ ^ ( ^ L 


-u [x 


f ) 


[x H , t) — \u [x H , t) | 


(x H ,i) 


u (x, 0) = 


~ £ ^X ~ 9 ( xL ' t ) =0 ’ 

( xH ’ t )-9( xH ,t)= 0, 

5 (0) (x) , 


(39) 


where u = u(x, t ) is the continuous solution vector, / = f(u) is a nonlinear flux function. The 
boundary conditions in (39) are constructed such that the senri-discrete energy only increases with 
respect to the imposed data and maintains the same form in the inviscid limit e — ► 0. 

A general semi-discretization of (39) suitable for energy or entropy analysis is 


A 

dt 


U 


+ 


+ 


’Zlgl^lg'D f (u) = Tlgl^lg'DS'Du 

u ( x L ) - e(Du) ( x L ) - g (x L )^ l LGL ^ LG V~ l e ( x L ) 

^ n(x ) j^(x )| u (x H ) - e(Du) (x 11 ) - g (x H )^j Ilgl^lgV^ 1 e (x H ) 


(40) 


where the initial data is u (x, 0) = g^(x) and f(u) denotes a numerical discretization of the inviscid 
flux (to be specified later). To simplify the notation in (40), the dependence on the time, t, has 
been omitted. 
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Two stability proofs for Burgers’ equation are now presented. The first proof uses conventional 
energy analysis and a canonical a-flux splitting technique, while the second proof uses entropy 
analysis. Both proofs assume that the solution is stored at the LG points and is interpolated to 
the LGL points to achieve a statement of stability. 


4.2.1 Energy analysis 

The flux /(«) in (39) can be written as the product of two functions, V(u) and W(u), (i.e. , 
f(u) = V(u)W(u)). Expressing the divergence term — in (39) as a combination of divergence 
and chain-rule forms, yields the following equation 


du df(u ) 

dt a dx 


+ (1 


a) 


V(u) 


dw{u) 

dx 


+ 


dV(u) 

dx 


W(u) 


o(^ iv> m) 

dx 


(41) 


Equation (41) is referred to as a-splitting and for a = 2/3 it is denoted a “canonical” splitting of the 
Burgers’ equation. Therefore, if we canonically split to the quadratic term in the semi-discretized 
Burgers’ equation (40) and apply the Sta-Grd-Alg to all spatial terms, we get the following staggered 
grid operator 


Au 

dt u 


+ 


g Zlgl^lg (P [n] + [n] V) u = 1 LG l^lg T>[e]Vu 

/ ^(x )+j^(.r )| u __ e (£> u ) (a ;£) _ g (x L )^j 1 L GL^LG e ( X L ) 

^ M (x ) ju(x )| u _ e (£> u ) X H ) _ g X LGL ^ LG -p - 1 e ( X H ) , 


(42) 


with the initial data u(x, 0) = g^(x). 


Theorem 4.1. The semi-discrete solution u defined in equation (42) is bounded for all time for 
any diagonal norm SBP operator V, provided there exist interpolation operators Ilg^lgl and 
Ilgl^lg that, satisfy the constraint 


VIlgl^lg = Ilg^lglP, 

and provided the boundary data \g {x L )\ and | g {x H )\ are bounded. 

Proof. The energy method is used to prove Theorem 4.1. The proof begins by contracting equation 
(42) with the vector u t T ) to yield the expression 


(43) 


| lu'Vu + 5 U 1 (Q [n] + [fl] Q) u = u 1 V[e\ V u 


+ 


u [x )+jtt(x )| ^ (x^ — e (Du) (x L ) — g (x L )^j u ( x L ) 
u { x ) j u (- r )l u X H _ e (x H ) — g (x H )^j u ( x H ) . 


The constraint given in Theorem 4.1 (see also equation (B29)) is used in all spatial terms to relate 
the vector contraction u t T ) on x, to an equivalent vector contraction on x. 


17 



Further simplification of equation (43) by using the SBP property Q+Q T = B and the boundary 
data yields the expression 

| ^u T ?uj = — e(Vu) T Wu — | (x H ) 3 — u (x L ) 3 ^j 

+ -u{x L f + \u(x H )\u{x H f -\u{x L )\u{x L f) (44) 

+ u (x L ) g (; x L ) + u ( x H ) g ( x H ) . 

Provided that u (x L ) / 0 or u ( x H ) / 0, the identity 

yz = ~2 -j= a Z ) + I y2 + ^ 2 ’ “ >0> (45) 

is used to bound the terms u ( x L ) g ( x L ) + u ( x H ) 5 (x ff ). The final expression becomes 
I (G T PG) = -e(Pu) T PPu + (,x L ) 2 + (: x H f 


\{Vla H )u ( x H ) - -j= g {x H )^j - \[yf{a L )u ( x L ) - ^3 = 5 (x L )) 


+ 1 v- 


(* H ) 2 +(£- 


u 


(* L ) 2 , 


(46) 


0 < a L < | | it (ar L ) I , 0<a*<§|«(®*)|. 

Thus, ^ |u T Puj < — £('Du) t V'Du and establishes boundedness of the L 2 norm of u provided that 
the boundary data g (x L ) 2 and g ( x H ) 2 are well behaved. □ 

An expression equivalent to (46) is now derived using entropy analysis [25]. 


4.2.2 Entropy analysis 

An entropy-entropy flux pair, (£, F), a potential-potential flux pair, (cj>, ?/>), and the entropy variable, 
w, for Burgers’ equation are [25] 

(S,F) = ; (fail’) = ; u = w. (47) 


Note that the entropy is guaranteed to be convex for all u (i.e., d 2 S/du 2 = 1), and that the entropy 
is (chosen) to be equivalent to the energy used in the SBP analysis. 7 

Consider the entropy analysis of equation (40). Apply the Sta-Grd-Alg to construct the 
quadratic inviscid and linear viscous fluxes as well as the boundary penalties. As with the col- 
located approach [1], the entropy and energy analyses of the time, viscous, and SAT terms are 
equivalent when using the Sta-Grd-Alg, while differences appear in the analysis of the quadratic 
flux term. In the entropy analysis, the quadratic flux is discretized using a diagonal norm SBP 
operator: T> = V~ l Q. The Q operator is then rearranged into telescoping form by using the gener- 
alized SBP relation given in Lemma 2.1. The resulting expression for the quadratic term discretized 
using Sta-Grd-Alg is 


1 du 2 

2 l)x 


V^Afiu). 


7 Tlie entropy is not unique. 
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The resulting semi-discrete staggered grid operator is then 


+ Ilgl^lg'P x Af(u) = Tlgl^lg F> [e] V u 

— E (x L ) - e (Vu) ( x L ) - g (x L )^j 1 LGL ^ LG V e ( x L ) 

u i x ) j^(- T )l tf _ e (£> u ) _ g (x H )^j 1 LGL ^ LG V e (x H ) 

with the initial data u(x,0) = g(°)(x). What remains is to construct an entropy conserving flux 

f( u ) = f! s) - 

Theorem 3.3 guarantees the existence of an entropy conserving flux for the delta form opera- 
tor 'P~ 1 Af(u), provided there exists a two-point entropy flux relation and a diagonal norm SBP 
operator is used for V. Using the definition for the entropy variable u given in equation (47) and 
Tadmor’s integral relation [25] 



fs ( u k , ue) = J 9 (w(u k ) + f (w(u e ) - w(u k ))) d£, g{w(u)) = f(u ) (49) 

o 


yields a two-point entropy conservative flux 


1 


f s u k ) — g ( u £ + u i u k + U k ) i 


which satisfies the two-point relation 


1 


(u £ - u k ) f s (u e , u k ) = - (uf - u\) . 

The high-order accurate entropy conserving flux guaranteed by Theorem 3.3 is given by 

N i N i 


(50) 


(51) 


7f } = E E 2 qekf S {ue,u k ) = 2 E E Qek^ (uj + ueu k + u k ) , 1 < i < N — 1. (52) 


1 


k=i -\- 1 t=\ 


k=i + 1 £=1 


6 


Contracting the quadratic term Zl G l^lg'P~ 1 At in equation (48) with the discrete vector 


T 


Vu ) yields the telescoping condition 

u T Af (S) = l T [n]Af (S) = 1 t AF = (Fjv-Fi) 


with 


N i 


Fi = EE qek {ui + u k ) f s (u e u k ) - - (u\ + u k ) 

k=i -\- 1 £= 1 ^ 


1 <i < N - 1, 


(53) 


(54) 


- 1 o — 1 

f n = 


Collecting all terms in the entropy analysis of Burgers’ equation yields an (entropy) estimate that 
is identical to the energy estimate given in equation (46). 
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4.2.3 Relating the energy and entropy analyses 

The canonical splitting of the inviscid flux ( f{u ) = = V{u)W{u)) in Burgers’ equation 

satisfies the necessary conditions of equation (A9). As such, all terms can be rearranged into a 
single telescoping flux 


1 du 2 df(u ) 

2 dx 8x 


+ (1 


a) 


V{n) 


dw{u) 

dx 


+ 


dV(u) 

dx 


W{u) 


P^ 1 Af(u). 


The a — splitting of the conservative (equation (55a)) and chain rule (equation (55b)) forms of the 
flux are decomposed as follows: 


V{u) 

V{u) 


-UU 1 
2 

; W(u) = 1 

(55a) 

1 Ml 

2 

; W(u) = U 1 

(55b) 


which produce the equivalent fluxes (see equation (A10)) 


N i / 

E W I UkUk + UeUg 

9 


-p a 

J i 


f\ = 


k=i -\- 1 1=1 
N i 


V 


EW ( UkU( + U^Uk \ 

9 ) 


k=i + 1 1=1 


1 < i < N - 1, 

1 < i < N- 1. 


(56a) 

(56b) 


Combining the split Burgers’ fluxes (56) into equation (A10) yields the expression 


fi = fi+fi = 2 i^qek^(u 2 e +u e u k + ul), 

ktt+iti 6 (57) 

1 < i < N - 1, / 0 = -u\, f N = 
for a splitting parameter a = 2/3. 

.Remark. Equation (40) invokes a canonical splitting of Burgers’ equation to facilitate the 
energy analysis, and is valid for any diagonal norm SBP operator. Equation (53) relies on a two- 
point dyadic entropy conservative flux of the form given as a necessary condition in Theorem 32. 

Choosing u as the entropy variable in equation (47), enforces equality of the entropy and a-split 

—(S) 

fluxes. (Compare the fluxes f i in equations (57) and (52).) 

Remark. The compressible Navier-Stokes equations do not support a canonical decomposition 
based on the a-split flux technique. Thus, conventional nonlinear energy analysis is not applicable 
(to our knowledge). Nevertheless, the existence of a two-point entropy conservative flux satisfying 
equation (32) enables entropy analysis to be used for all diagonal norm SBP operators. 

Remark. The numerical approach is conservative. The inviscid and viscous fluxes are explicitly 
conserved on the LGL points. 
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5 The compressible Navier-Stokes equations in multiple dimen- 
sions on staggered grids 

5.1 Staggered grid in two dimensions 


Extension of a ID staggered operator to multiple dimensions can proceed in several ways. Figure 3 
shows two popular staggered data structures in two spatial dimensions. Both approaches store the 
solution at the tensor product LG points (i.e. , the so called Gauss points) of order p (blue crosses). 

The fully-staggered approach moved the data via 3D tensor product interpolations from LG 
(supporting a polynomial of order p) to LGL (supporting a polynomial of order p+ 1) points (black 
circles). The discrete operators reported in references 1,2,13,30, are used on the LGL points to 
construct the spatial residual. The temporal updates needed on the LG points are obtained by 
restricting the LGL residuals back to the LG points. Extension to general curvilinear coordinates 
follows immediately on the LGL points [1,13]. 

The semi-staggered approach moves the data via ID interpolations from the LG to LGL points 
(black circles and green triangles). The inviscid terms are constructed via three ID operations on 
the semi-staggered LG-LGL points. The viscous terms are most easily formed by using a fully- 
staggered approach. The semi-staggered operator has the advantage of not requiring corner data 
in the inviscid operators. However, its extension to curvilinear coordinates is not straight forward 
because of ambiguities in the geometric conservation law (GCL) terms. For this reason, the fully- 
staggered approach is used exclusively herein. 
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(a) Fully-staggered. 


(b) Semi-staggered. 


Figure 3. Fully- and semi-staggered 2D tensor product elements. 
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5.2 Tensor operators in three dimensions 

Consider a single tensor product element and an entropy stable spatially discontinuous collocation 
(SSDC) discretization with M = p+ 1 LG solution points in each coordinate direction 8 [1,2,30,31]; 
the following element-wise matrices will be used: 

P = (Vm 0 Pm 0 Pm 0 h^j , 


ZlG^LGL = {PlG^LGl) XiX2X 3 = {PlG^LGL 0PlG^LGL 0PlG^LGL 0 h) , 
ZlGL^LG = {PlGL^Lg) X1 X 2 X 3 = {Plgl^lg 0Plgl^lg 0Plgl^lg 0 h ) , 
V Xl = ( T>n 0 In 0 In 0 h) , ■ ■ • Px 3 = ( In 0 In 0 Pn 0 h) , 

Pxi = {Pn 0 In 0 In 0 h) , ■ ■ ■ Px 3 = {In 0 In 0 Pn 0 h) , 

P x ix 2 = (Pn 0 Pn 0 In 0 h) , Px 2 x 3 = (In 0 Pn 0 Pn 0 h) , 

P = Px ix 2 x 3 = (Pn 0 Pn 0 Pn 0 h) , 

B X1 = (Bn 0 In 0 In 0 h) , ■ - = (/at 0 In 0 Bn 0 h ) , 

An = (Aat 0 In 0 In 0 h) , ■■■ A X3 = (Jjv 0 In 0 An 0 h) , 


where Pm is the norm of the LG points, while V at, Pn , Aat, and Bn are the ID SBP operators [31] 
defined on the LGL points, and In is the identity matrix of dimension N. I 5 denotes the identity 
matrix of dimension five . 9 The subscripts in (58) indicate the coordinate directions to which the 
operators apply (e.g., T> X i is the differentiation matrix in the x\ direction). The symbol <g> represents 
the Kronecker product. When applying these operators to the scalar entropy equation in space at 
the LG points, a hat is used to differentiate the scalar operator from the full vector operator. For 
example, 

P = (p m 0 Pm 0 Pm ) ; P = (Pm 0 Pm 0 Pm) ■ (59) 

The vector of conservative variables of each element is ordered as 

q= (9 (®(1)(1)(1)) T ,<7 0c(i)(i)(2)) T , ■ ■ ■ ,q (x(M)(M)(m)) T ) = (?(1),<Z(2), • • • .9(M3)) > ( 60 ) 

s Recall from Section 2.2.2 that with the staggered algorithm the number of LG and LGL points in ID is denoted 
by M and N, respectively. 

9 The 3D compressible Navier-Stokes equations form a system of five nonlinear partial differential equations. 
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where the subscripts denote the ordering of the solution points in the coordinate directions. Assume 
an equivalent definition and order for the entropy variables w; and an analogous definition for the 
variables at the LGL points. 10 

The Zlg^lgl operator transfers data from x to x, while the Tlgl-*lg transfers data from x 
to x. Define 


W = ; Cij = iLG^LGLCij ; [cij] = Diag [l LG -*LGLCij]- (61) 


Using these definitions and the SBP operators (58), system (18) is discretized locally on an 
isolated element as [13, 31] 


dq 

dt 


+ Zlgl^lg 


r, 1 a ,, f i - V x ; f 


Xi 


n(VY 

Xi M 


= Zlgl^lg V Xi ] 


(62) 


where Einstein notation is used to express the coordinate directions. The penalty interface terms 
g( Int ^ with i = 1,2,3 are used to connect neighboring elements (see Section 6). 

An entropy conservative reconstruction is used for the inviscid fluxes fj in equation (62) based 
on the expression given in theorem 3.3 and the two-point entropy conservative flux, f s (u£,Uk), 
of Ismail and Roe [29]. The interpolated entropy variables w = Tlg^lgl w are used to build 
the fluxes on the LGL points. The viscous fluxes are also computed using interpolated entropy 
variables and the operators T> Xi , i = 1,2,3, defined in (58). The viscous coefficient matrices c^- are 
again formed using interpolated data on the LGL points. 


Remark 5.1. The interpolations from and to the LG points are carried out in computational space 
by using an efficient tens or- product algorithm that requires only the knowledge of the ID Tlg^lgl 
and Ilgl^lg operators. The extension to general curvilinear coordinates follows immediately on 
the LGL points [1, 2j. 


6 Entropy stable interface coupling 


Consider two cubic tensor product elements by extending equation (62) to two adjoining elements. 
Without loss of generality assume that all their faces are orthogonal to the three coordinate direc- 
tions and are not boundary faces, i.e. , they are not part of the boundary surface dQ. The resulting 
expressions become [13, 31] 


dqe 

~dt + -^ LGL ^ LGi 


'P 'xi ,t ^ x i / D Xi / [cud 


'■T') — 1 (-/"?!<£) 

= l L GL^LG(.l J Xii £ g,:y 


©M-^w^^g^’ 0 , 


dq r 

dt 


+ ZlGL^LGt [D xi,r A Xi,r f i,r D Xi ^ r [Cjj,r] ®j,r_ ~ -^LGL^LGr-D ' Xi , r § 


-1 

i,r ’ 


Q ir - v Xi w i = p x .] r g£? t),e t 


(63a) 

(63b) 

(63c) 

(63d) 


where the subscripts l and r denote the “left” and “right” elements, ©n and ©i <r are the vectors 
of the gradient of the entropy variables on the left and right elements in the i direction, whereas 


1() With the staggered algorithm, the number of LGL points in ID is denoted by N (see Section 2.2.2). Thus, N 3 is 
the number of LGL points in a 3D tensor-product element. 
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g anc ^ §.■ are the penalty interface terms on the conservative variable and the gradient 
of the entropy variable, respectively [31]. As indicated in (61), the matrices [c^] are block diagonal 
matrices with IV 3 five-by-five blocks corresponding to the viscous coefficients of each LGL point. 11 
Note that (63) is obtained by using j ' = c l3 w Xj = Cij Qj. 

To obtain an equation for the entropy of the system, we follow the entropy stability analysis 
presented in [1,31]. Therefore, multiplying the two discrete equations in the left element by w JVi 
and V/, respectively, and the two discrete equations in the right element by wj V r and 

([%,r]©j.r) T T’r, respectively, the expression for the time derivative of the entropy function S in 
each element is 


d_ 

dt,~ 


^1 T V, Si+ 2 


. c ij,l\ ®j,l 


+ 1 ( 'Px 2 X3,l B Xl ,l*l,l + 'Pxixzt B X2 , 1*2,1 + 'Px\x 2 ,1 B X3 ,l*3,l 


W l (V X2X3 ,1 B X\ ,1 T V x ,X3,l ^X2,l \p2 j,l\®j,l T V x\X2,l ^X3,l [Gy,/] ® ji,Z ) 


+ w 7 (v ,J lnt )’ q + v ,j lnt ^ q + v 

' W Z \ 1 X2X3 ,1 &1 l W r X \X3,l &2 t l I 7 X1X2 ,1 &3 t l J 


\T -rt (Jni) ,© 


+ ([c W \&u) ’ V X2 X3,l Z\7 h ” + V X1 X3 ,1 g£ nt) ’° + ([c m ]&,,l) T V X1X2 ,l g£ nt) ’ 0 > 

(64a) 


4l T Pr Sr + 2 

at 


© j,r 


~\~ 1 ( 'Px 2 x 3 ,r Bxi,r*l,r ^Px\x 3 ,r Bx 2 ,r*2,r 'Pxix 2 ,r Bx 3 ,r*3,r 


— W r (J^x 2 X 3 ,r B X i ,r [cij,r\®j,r ^Pxix 3 ,r B X2 ,r \^2j,r\®j,r T^x\x 2 ,r B x ^,r [^3 j,r\®j,r) 

i W T (V (/nt),# v ( Int),q v (Int),q\ 

i w r [' x 2 x 3 ,r &l 5 r ' ’ o2,r ' ' #i# 2 ,r& 3 jr J 

+ ([cii,r]©j,r) T ^ X2X3,r + ([%r] ©,,r-) T gg^’® + ([%>] ®,>) T V XlX2 , r gg? 0 ’®, 

(64b) 

where the vectors 1 and 1 represent a vectors with M 3 and IV 3 elements, respectively (i.e., 1 = 

(1, 1; > 1)m 3 )- 


To simplify the notation, assume that the interface between the two tensor product cells lies at 
x\ = 0. We also assume that all the points that lie on the other faces of the two cubes are treated 
in an entropy stable fashion; their contribution can then be neglected without loss of generality. 
Then, for our analysis, we can just focus on a pair of LGL interface nodes at x\ = 0. We then 
introduce the operators and eW, which “extract” from the cell- wise solution and flux vectors 
only the variables associated to these LGL points (i.e., the node at the left, (— ), and at the right, 
(+), of the interface). 12 Therefore, equations (64) reduce to 


d t l T ViSi + l T V x 2X3 ,iF hl + 2 


. 2 

= ™JV X2 X3,1 +wJr X2X3 ,ig[ I J lt),q 

' Pi 

(65a) 


+ {[ci j A@ j ,i) T v X2X 3 ,zgS nt) ’ e , 


11 In the staggered algorithm framework, N 3 is the number of LGL points in a 3D tensor-product element. 
12 These two vectors are zero at all points except at the “(— ), (+)” interface points. 
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2 

= -W Jv x 2 x 3 ,r [cij jr ]©j !r e (+) +wjv x 2 x 3 ,rgir t),q 

Vr 

+ ([cij,r]©j>) 'Px2X 3 ,rSi j r ■ 

(65b) 

The interface penalty terms are constructed as a combination of a local discontinuous Galerkin-type 
(LDG-type) approach and an interior penalty (IP) technique [13]: 



The LDG penalty terms involve the coefficients ^(1 ±q) and act only in the normal direction 
to the face. The IP terms involve the block diagonal parameter matrix, [L] = Diag [L ] , with N 3 
five-by-five blocks, L , which are left unspecified for the moment. 13 

Herein, the solution between adjoining elements is allowed to be discontinuous. An inviscid 
interface flux that preserves the entropy consistency of the interior high-order accurate spatial 
operators [1] on either side of the interface f lySS ' 1 (q^~\q^) is constructed as 

f(SS) ^(-) J ,(+)) = f{SC) + A («,(+> - «;H) , (67) 

where f^ sc (q^~\ is the entropy conservative inviscid interface flux of any order [1,2,19,30,31] 
—(S) 

(i.e., / ). This flux is constructed as in Theorem 3.3 by using the two-point entropy conservative 

flux, f s (u£,Uk), of Ismail and Roe [29]. A is a negative semi-definite interface matrix with zero 
or negative eigenvalues. The superscripts (— ) and (+) denote the collocated values on the left 
and right side of the interface, respectively. The entropy stable flux f^ ss ^ (q^~\ q^) is more 
dissipative than the entropy conservative inviscid flux f^ sc ^ {q^~\q^), as can be easily verified 
by contracting f^ ss ^ (q^~\q^) against the entropy variables [1]. Note that in reference 15, grid 
interfaces for entropy stable finite difference schemes are studied and interface fluxes similar to (67) 
are proposed. 

13 In the staggered algorithm framework, N 3 is the number of LGL points in a 3D tensor-product element. 


-l'VrS r -l'V X 2 X 3 ,r) 
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Substituting expressions (66) and (67) in (65) and summing all the contributions of the two 
elements results in 


^l T ViSi + ^-l T V r S r + 2 


dt 


dt 


n, 

2 


2 


+ 

Vi 

\J [ c ij,r\ ®j,r 

Vr_ 


= T (/) + T (y) , (68) 


where Y^ and Y (l ) are the inviscid and the viscous interface terms. At the two interface nodes, 
these terms are 


= ^ w (+) — w ( f( ss ' > — il>( ^ — w/ ^ A — k/ 

(69) 

Y (V •> = L (u; (+) — . (70) 

Clearly, the interface contributions are dissipative if both the five-by-five matrices A and L are 
negative semi-definite. The matrix A can be constructed using different approaches, e.g., using an 
upwind operator that dissipates each characteristic wave based on the magnitude of its eigenvalue: 


f ssc (? (_) ,9 (+) ) = f (SC) (q ( ~\q {+) ) + i/2Y|A| y T (w {+) - , 


d_ 

8q 

dq 


/(?) 


dw 


= yxy T , 
yy T , 


(71) 


where A and y are the diagonal matrix of the eigenvalues and the matrix of the eigenvectors, 
respectively. Note that the relation dq/dw = TT T is achieved by an appropriate scaling of 
the rotation eigenvectors. Reference 1 constructs the matrix dq/dw using the scaled eigenvectors 
introduced by Merriam [32], The imposed artificial viscosity satisfies the senri-discrete second law 
of thermodynamics. 

We are then left with the viscous interface term, Y ! ' j . The parameter values a = 0 and 
a = ±1 yield a symmetric LDG and a “flip-flop” narrow stencil (nearest neighbor) LDG penalty, 
respectively [2,31]. An LDG value of a = 0 produces a global discrete operator that has a neutrally 
damped spurious eigennrode. Herein, when a = 0, this mode is damped using the IP dissipation. 
For a Reynolds number that approaches oo, we would like the five- by- five matrix L to go to zero so 
that only the inviscid entropy stable penalty contributions in (66) are recovered. To achieve that, 
the matrix L is constructed as 


, x ~ (-) + ? (+) 


2 (V. 


/ 3 ( Int) > 0 , 


(72) 


where c^ 1 ' ) and are the positive semi-definite viscous coefficient matrices at the left and right 
side of the interface in the normal direction. The coefficient ) can be used to modify the 
strength of the IP penalty term, although excessively large values of /3( In d reduce the maximum 
stable time step. Reference [2, 31] selects to be the maximum value for which the explicit 

stability constraint remains unaffected. The factor (^ > x 1 )( 1 )( 1 ) i n the denominator represents the 
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normal local grid spacing, which is incorporated in the diagonal SBP operator TT 14 This term is 
introduced to get the correct dimension, and, as for the standard IP finite element approach, it 
increases the strength of L with increased resolution. 


7 Discussion: A theoretical cost analysis 


A p - th order 3D tensor-product LGL element has (p+ l) 3 collocation points. The total cost of 
the fully discrete, explicit-in-tinre operator, predominantly results from the inviscid and viscous 
differentiation operators. Thus, the operation count for a conventional non-staggered element [1,2] 
scales as (p + l) 4 . By a similar accounting, the operation count for the semi-staggered operator 
scales as (p + 1 )(p + 2) 2 and (p + 2) 4 for the fully-staggered operator. The work ratio between the 
conventional algorithms and the two staggered approaches is given by 


p + 2 
p+ 1 


S 


= Wgt 


(73) 


with s = 2 (semi-staggered), and s = 4 (fully-staggered), respectively (assuming that the explicit 
time stepping stability constraint is equivalent for all schemes). 

The results shown in Section 8 demonstrate that the fully-staggered SSDC approach is signifi- 
cantly more accurate than the conventional LGL SSDC operator for the same solution polynomial 
and grid density. To determine which approach is more efficient, assume that two simulations 
are performed on the same grid, and that r e is the ratio of errors between the conventional and 
fully-staggered operators. A comparable accuracy can be achieved by increasing the number of 
elements used in the conventional approach. Thus, assume that N^ onv and Ng tag are the number 
of elements required by each approach to achieve the same error. The work ratio ^Conv/^Stag f° r 
ID elements of polynomial order p. required to overcome the differential in accuracy, is therefore 

r e = [N Conv/^Stag ) • N° w assume that the error is homogeneously distributed in D spatial 

dimensions, and that a decrease in element size produces a proportional change in the explicit 
temporal stability constraint. The resulting expression for the cost ratio is given by 

r>+l Nn 

(r e )L+T = -^onv = Whr (74) 

iV Stag 

Figure 4 shows the ratio of cost between the two approaches as a function of the polynomial 
order p. Five error ratios are assumed: r e = 2^,0 < (3 < 4. The first comparison (r e = 1) assumes 
the two approaches to be of comparable accuracy. The last r e = 16 assumes that the staggered 
operator is 16 times more accurate. This simple estimate suggests that the fully staggered approach 
is cost effective, despite the increased operation count. 


8 Numerical results: Accuracy and robustness 

8.1 Isentropic Euler vortex propagation 

The test case considered in this section is the (inviscid) propagation of a vortex for which an exact 
solution is known. This is an excellent test problem for verifying the accuracy and functionality of 

14 Herein, with i = 1 appears in the definition of L because the interface between the two elements is 

orthogonal to the xi direction [13,31]. 
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Figure 4. Cost comparison of conventional LGL [1,2] and fully-staggered SSDC operators. 


the inviscid components of a compressible Navier-Stokes solver. It is fully described by 

f(xi,x 2 ,x 3 ,t) = 1 - (x! - xgo - Uoocos(a) t ) 2 + (x 2 - x 2 ,o - t/ooSin(a) t ) 2 

T(xi,x 2 ,x 3 ,t) = 1 - elM^ 0 ^exp(f(xi,x 2 ,x 3 ,t)) 

1 

p{xi,X 2 ,X 3 ,t) = Tt-1, 

/ j \ tt r \ X2—X2,0—Uoos'm(a)t ( f(x\,X2,X3,t ) ' 

rii(xi,x 2 ,x 3 ,f) = C/oo cos(a) — e t , —exp 1 


u 2 {xi,X2,x 3 ,t) = Uoo sin(a) - J ' 1,u 2 ^°° COb(aK exp ( ; 

U 3 (x 1 ,X2,X 3 ,t) = 0 , 


(75) 


where Goo, and (xgo, x 2 ) o, ® 3 ,o) are the module of the freestream velocity, the Mach number, 
and the coordinates of the vortex center, respectively. In this study, the values = ATocCoc, 

e v = 5.0, Mqo = 0.5, 7 = 1.4, and a = 45° are used; and the domain is described by 


me (-5,5), x 2 e ( 5, 5) , x 3 e(-i,i), (xi ) 0 ,x 2 ,o,m,o) = (0,0,0), t > 0. 

The boundary conditions are prescribed by penalizing the numerical solution against the exact 
solution by using an SAT approach [33], whereas the interface coupling between two adjoining 
elements is imposed using the entropy stable treatment proposed in references 1,13. 

The accuracy of the following entropy stable schemes is investigated using uniform Cartesian 
and unstructured nonuniform grids: 

• Fully-staggered SSDC algorithm with plg = 1,2, 3, 4, 5, 10 and plgl = 2, 3, 4, 5, 6 , 11. 

• Conventional LGL SSDC algorithm [1,13] with plgl = 1,2, 3, 4, 5, 10. 


8.1.1 Uniform Cartesian grid 

Different grid resolutions are examined, and the vortex is halfway out of the domain when the 
error measure is evaluated. This measures the effect of the penalty boundary conditions and the 


28 


interior scheme. Tables 1-6 show the convergence study for a sequence of maximum nine nested 
grids. 15 The number of cells in each coordinate direction is indicated in the first column of these 
tables (i.e. , “Resolution”). The L 2 norm of the error decay asymptotes towards the designed rate 
in each case (i.e., second order, third order, fourth order, fifth order, sixth order, and eleventh 
order, respectively), for both fully staggered and the conventional SSDC approaches. We highlight 
that the fully-staggered approach is significantly more accurate than the conventional algorithm for 
the same grid resolution. However, simple counting arguments based on inviscid flux evaluations, 
indicate that the cost of the staggered algorithm for a solution polynomial order p is comparable 
to that of an LGL operator [1,2], with a solution polynomial order of (p+ 1). 


Table 1 Error convergence is shown for the fully-staggered plg = 1> Plgl = 2 and conventional 
LGL [1, 2] plgl = 1 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, 

Plg = 1) Plgl 

= 2 

Conventional, plgl = 1 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.42e-02 

- 

6.28e-02 

- 

9.78e-03 

- 

3.70e-02 

- 

4x4x2 

8.73e-03 

0.70 

6.51e-02 

-0.05 

7.75e-03 

0.33 

5.11e-02 

-0.46 

8x8x2 

3.86e-03 

1.18 

3.75e-02 

0.80 

7.41e-03 

0.07 

6.33e-02 

-0.31 

16 x 16 x 2 

1.49e-03 

1.37 

1.52e-02 

1.30 

3.95e-03 

0.91 

4.08e-02 

0.63 

32 x 32 x 2 

5.35e-04 

1.47 

6.07e-03 

1.32 

1.66e-03 

1.25 

1.88e-02 

1.12 

64 x 64 x 2 

1.72e-04 

1.64 

1.95e-03 

1.64 

6.23e-04 

1.42 

7.38e-03 

1.35 

128 x 128 x 2 

4.76e-05 

1.85 

5.30e-04 

1.88 

2.12e-04 

1.56 

2.49e-03 

1.57 

256 x 256 x 2 

1.14e-05 

2.06 

1.27e-04 

2.06 

6.10e-05 

1.79 

6.96e-04 

1.84 

512 x 512 x 2 

2.72e-06 

2.07 

3.22e-05 

1.98 

1.46e-05 

2.07 

1.65e-04 

2.08 


Table 2 Error convergence is shown for the fully-staggered plg = 2, plgl = 3 and conventional 
LGL [1, 2] plgl = 2 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, 

Plg = 2, plgl 

= 3 

Conventional, plgl = 2 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.31e-02 

- 

9.28e-02 

- 

6.57e-03 

- 

6.44e-02 

- 

4x4x2 

3.52e-03 

1.90 

2.20e-02 

2.07 

6.14e-03 

0.10 

7.24e-02 

-0.17 

8x8x2 

8.83e-04 

2.00 

8.38e-03 

1.40 

2.66e-03 

1.21 

3.11e-02 

1.22 

16 x 16 x 2 

1.79e-04 

2.30 

2.78e-03 

1.59 

6.91e-04 

1.94 

1.13e-02 

1.46 

32 x 32 x 2 

3.10e-05 

2.53 

5.54e-04 

2.33 

2.01e-04 

1.78 

5.24e-03 

1.11 

64 x 64 x 2 

5.30e-06 

2.55 

8.86e-05 

2.64 

2.66e-05 

2.92 

7.15e-04 

2.87 

128 x 128 x 2 

7.80e-07 

2.77 

1.32e-05 

2.74 

4.35e-06 

2.61 

1.15e-04 

2.64 

256 x 256 x 2 

1.14e-07 

2.78 

2.19e-06 

2.60 

6.57e-07 

2.73 

1.75e-05 

2.71 

512 x 512 x 2 

1.42e-08 

3.00 

2.67e-07 

3.03 

8.57e-08 

2.94 

2.60e-06 

2.75 


15 The number of grid cells is doubled every time in each coordinate direction. 
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Table 3 Error convergence is shown for the fully-staggered plg = 3, plgl = 4 and conventional 
LGL [1,2] plgl = 3 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, 

Plg = 3, plgl 

= 4 

Conventional, plgl = 3 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

5.63e-03 

- 

3.30e-02 

- 

6.94e-03 

- 

8.21e-02 

- 

4x4x2 

1.77e-03 

1.67 

2.12e-02 

0.64 

3.69e-03 

0.91 

5.18e-02 

0.66 

8x8x2 

2.02e-04 

3.13 

3.60e-03 

2.56 

6.79e-04 

2.44 

1.12e-02 

2.21 

16 x 16 x 2 

1.71e-05 

3.56 

3.22e-04 

3.48 

7.30e-05 

3.22 

1.76e-03 

2.67 

32 x 32 x 2 

1.36e-06 

3.66 

2.76e-05 

3.54 

7.10e-06 

3.36 

1.60e-04 

3.46 

64 x 64 x 2 

1.01e-07 

3.75 

1.84e-06 

3.91 

7.07e-07 

3.33 

1.94e-05 

3.05 

128 x 128 x 2 

6.31e-09 

4.00 

1.15e-07 

4.00 

5.41e-08 

3.71 

1.80e-06 

3.43 

256 x 256 x 2 

3.51e-10 

4.17 

6.41e-09 

4.16 

3.69e-09 

3.88 

1.45e-07 

3.63 

512 x 512 x 2 

2.12e-ll 

4.05 

4.48e-10 

3.84 

2.40e-10 

3.94 

1.12e-08 

3.69 


Table 4 Error convergence is shown for the fully-staggered plg = 4, plgl = 5 and conventional 
LGL [1, 2] plgl = 4 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, 

Plg = 4, plgl 

= 5 

Conventional, plgl = 

4 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

4.14e-03 

- 

3.61e-02 

- 

5.76e-03 

- 

5.92e-02 

- 

4x4x2 

6.04e-04 

2.78 

8.41e-03 

2.10 

1.41e-03 

2.03 

2.05e-02 

1.53 

8x8x2 

3.17e-05 

4.25 

3.37e-04 

4.64 

1.17e-04 

3.59 

1.62e-03 

3.66 

16 x 16 x 2 

1.68e-06 

4.23 

3.45e-05 

3.29 

7.24e-06 

4.02 

1.79e-04 

3.18 

32 x 32 x 2 

6.40e-08 

4.72 

1.35e-06 

4.68 

3.91e-07 

4.21 

1.20e-05 

3.90 

64 x 64 x 2 

2.26e-09 

4.82 

5.39e-08 

4.65 

1.76e-08 

4.47 

5.57e-07 

4.43 

128 x 128 x 2 

7.36e-ll 

4.94 

1.99e-09 

4.76 

7.16e-10 

4.62 

2.41e-08 

4.53 

256 x 256 x 2 

- 

- 

- 

- 

2.19e-ll 

5.03 

8.37e-10 

4.85 


Table 5 Error convergence is shown for the fully-staggered plg = 5, plgl = 6 and conventional 
LGL [1, 2] plgl = 5 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, 

Plg = 5, plgl 

= 6 

Conventional, plgl = 5 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

2.73e-03 

- 

3.49e-02 

- 

3.05e-03 

- 

5.75e-02 

- 

4x4x2 

1.81e-04 

3.91 

1.34e-03 

4.71 

5.25e-04 

2.54 

5.62e-03 

3.36 

8x8x2 

8.68e-06 

4.38 

1.98e-04 

2.75 

2.60e-05 

4.33 

7.85e-04 

2.84 

16 x 16 x 2 

1.25e-07 

6.12 

2.76e-06 

6.16 

7.52e-07 

5.11 

2.46e-05 

5.00 

32 x 32 x 2 

2.63e-09 

5.57 

5.60e-08 

5.63 

1.90e-08 

5.31 

6.30e-07 

5.29 

64 x 64 x 2 

4.18e-ll 

5.98 

1.00e-09 

5.80 

4.00e-10 

5.57 

1.81e-08 

5.12 

128 x 128 x 2 

- 

- 

- 

- 

7.67e-12 

5.71 

3.71e-10 

5.61 
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Table 6 Error convergence is shown for the fully-staggered plg = 10, plgl = H and conventional 
LGL [1,2] plgl = 10 SSDC algorithms for the isentropic Euler vortex propagation on uniform 
Cartesian grids. 



Fully-staggered, p LG = 10, plgl = H 

Conventional, plgl = 

10 

Resolution 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.23e-04 

- 

1.36e-03 

- 

1.33e-04 

- 

2.41e-03 

- 

4x4x2 

4.35e-07 

8.15 

5.75e-06 

7.88 

1.83e-06 

6.18 

3.14e-05 

6.26 

8x8x2 

6.17e-10 

9.46 

1.79e-08 

8.33 

3.30e-09 

9.12 

l.lle-07 

8.15 

16 x 16 x 2 

2.07e-12 

8.22 

1.57e-ll 

10.15 

1.27e-12 

11.35 

5.28e-ll 

11.03 
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8.1.2 Unstructured nonuniform grids 


The goal of this study is to investigate whether the fully-staggered algorithm is effectively (p+ 1)- 
order accurate for more realistic meshes (p is the order of the polynomial supported by the LG 
points). As for the case of uniform Cartesian grids, different grid resolutions are examined, and 
the vortex is halfway out of the domain when the error measure is evaluated. A nested family of 
grids is constructed by replicating N times the “unstructured grid kernel” shown in Figure 5. 



Figure 5. “Unstructured grid kernel” used to construct a sequence of nested grids for the inviscid 
vortex test case. 

Tables 7-12 show the convergence study for a sequence of twelve nested grids. The number of 
“unstructured grid kernels” in each coordinate direction is indicated in the first column of these 
tables. For the fully-staggered SSDC algorithm, the L 2 norm of the error decay asymptotes towards 
the designed rate (i.e., (p + 1)) and, in each case, is more accurate than the conventional LGL SSDC 
algorithm for the same grid resolution. The conventional path converges instead to p. 
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Table 7 Error convergence is shown for the fully-staggered plg = 1> Plgl = 2 and conventional 
LGL [1,2] plgl = 1 SSDC algorithms for the isentropic Euler vortex propagation on highly nonuni- 
form grids. 



Fully-staggered, 

Plg = 1, Plgl 

= 2 

Conventional, plgl = 

1 

Resolution 

L ? error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

1.19e-02 

- 

1.34e-01 

- 

8.22e-03 

- 

8.77e-02 

- 

2x2x2 

6.31e-03 

0.91 

8.16e-02 

0.72 

1.07e-02 

-0.38 

2.20e-01 

-1.33 

3x3x2 

5.12e-03 

0.51 

7.39e-02 

0.25 

1.02e-02 

0.11 

1.37e-01 

1.18 

4x4x2 

3.57e-03 

1.25 

5.52e-02 

1.01 

1.07e-02 

-0.14 

1.94e-01 

-1.21 

5x5x2 

2.56e-03 

1.49 

4.94e-02 

0.50 

1.02e-02 

0.21 

2.02e-01 

-0.18 

6x6x2 

1.87e-03 

1.73 

2.95e-02 

2.83 

8.82e-03 

0.78 

1.70e-01 

0.96 

7x7x2 

1.52e-03 

1.36 

2.91e-02 

0.08 

7.36e-03 

1.18 

1.46e-01 

0.98 

8x8x2 

1.26e-03 

1.39 

2.36e-02 

1.58 

6.36e-03 

1.09 

1.14e-01 

1.81 

9x9x2 

1.04e-03 

1.61 

2.27e-02 

0.33 

5.82e-03 

0.75 

1.08e-01 

0.47 

10 x 10 x 2 

8.76e-04 

1.65 

2.06e-02 

0.93 

5.43e-03 

0.66 

1.06e-01 

0.20 

11 x 11 x 2 

7.52e-04 

1.61 

1.69e-02 

2.06 

5.06e-03 

0.75 

1.02e-01 

0.38 

12 x 12 x 2 

6.54e-04 

1.59 

1.60e-02 

0.63 

4.70e-03 

0.84 

1.16e-01 

-1.43 

13 x 13 x 2 

5.75e-04 

1.61 

1.62e-02 

-0.16 

4.39e-03 

0.87 

1.13e-01 

0.29 

14 x 14 x 2 

5.09e-04 

1.67 

1.45e-02 

1.54 

4.11e-03 

0.87 

1.05e-01 

1.06 

15 x 15 x 2 

4.53e-04 

1.66 

1.27e-02 

1.83 

3.87e-03 

0.88 

9.40e-02 

1.56 

16 x 16 x 2 

4.07e-04 

1.66 

1.17e-02 

1.33 

3.65e-03 

0.90 

8.30e-02 

1.92 


Table 8 Error convergence is shown for the fully-staggered plg = 2, Plgl = 3 and conventional 
LGL [1,2] plgl = 2 SSDC algorithms for the isentropic Euler vortex propagation on highly nonuni- 
form grids. 



Fully-staggered, 

Plg = 2, plgl 

= 3 

Conventional, plgl = 

2 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

5.25e-03 

- 

3.96e-02 

- 

9.21e-03 

- 

1.47e-01 

- 

2x2x2 

1.62e-03 

1.70 

1.88e-02 

1.07 

3.95e-03 

1.22 

7.74e-02 

0.92 

3x3x2 

1.38e-03 

0.39 

2.28e-02 

-0.47 

3.14e-03 

0.56 

6.92e-02 

0.28 

4x4x2 

7.71e-04 

2.03 

2.18e-02 

0.15 

2.53e-03 

0.76 

7.05e-02 

-0.07 

5x5x2 

3.21e-04 

3.92 

7.48e-03 

4.80 

1.65e-03 

1.90 

5.91e-02 

0.79 

6x6x2 

2.04e-04 

2.49 

4.28e-03 

3.07 

1.14e-03 

2.06 

3.37e-02 

3.08 

7x7x2 

1.47e-04 

2.11 

3.35e-03 

1.58 

8.75e-04 

1.70 

2.81e-02 

1.19 

8x8x2 

1.08e-04 

2.34 

2.87e-03 

1.15 

6.86e-04 

1.82 

2.14e-02 

2.02 

9x9x2 

8.10e-05 

2.43 

2.14e-03 

2.49 

5.37e-04 

2.08 

1.75e-02 

1.73 

10 x 10 x 2 

6.18e-05 

2.57 

1.56e-03 

3.01 

4.26e-04 

2.19 

1.59e-02 

0.90 

11 x 11 x 2 

4.76e-05 

2.75 

1.41e-03 

1.06 

3.50e-04 

2.08 

1.21e-02 

2.82 

12 x 12 x 2 

3.71e-05 

2.86 

1.27e-03 

1.20 

2.95e-04 

1.95 

9.28e-03 

3.08 

13 x 13 x 2 

2.95e-05 

2.85 

8.69e-04 

4.75 

2.54e-04 

1.88 

7.92e-03 

1.97 

14 x 14 x 2 

2.39e-05 

2.85 

6.44e-04 

4.05 

2.21e-04 

1.84 

7.48e-03 

0.78 

15 x 15 x 2 

1.97e-05 

2.82 

5.82e-04 

1.45 

1.95e-04 

1.82 

7.05e-03 

0.85 

16 x 16 x 2 

1.64e-05 

2.80 

4.59e-04 

3.67 

1.74e-04 

1.81 

6.48e-03 

1.31 
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Table 9 Error convergence is shown for the fully-staggered plg = 3, plgl = 4 and conventional 
LGL [1,2] plgl = 3 SSDC algorithms for the isentropic Euler vortex propagation on highly nonuni- 
forrn grids. 



Fully-staggered, 

Plg = 3, plgl = 4 

Conventional, plgl = 

3 

Resolution 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

2.33e-03 

- 

3.57e-02 

- 

4.76e-03 

- 

5.64e-02 

- 

2x2x2 

4.73e-04 

2.30 

6.71e-03 

2.41 

1.02e-03 

2.22 

2.04e-02 

1.47 

3x3x2 

2.39e-04 

1.69 

7.59e-03 

-0.31 

9.00e-04 

0.32 

2.19e-02 

-0.18 

4x4x2 

8.94e-05 

3.42 

2.14e-03 

4.40 

4.70e-04 

2.26 

1.87e-02 

0.54 

5x5x2 

4.59e-05 

2.99 

1.27e-03 

2.34 

2.79e-04 

2.34 

1.05e-02 

2.59 

6x6x2 

2.27e-05 

3.86 

6.72e-04 

3.49 

1.67e-04 

2.80 

7.92e-03 

1.55 

7x7x2 

1.32e-05 

3.51 

4.66e-04 

2.37 

9.27e-05 

3.83 

4.74e-03 

3.33 

8x8x2 

8.58e-06 

3.24 

4.55e-04 

0.18 

5.82e-05 

3.49 

3.00e-03 

3.43 

9x9x2 

5.76e-06 

3.39 

3.41e-04 

2.44 

4.05e-05 

3.07 

1.82e-03 

4.23 

10 x 10 x 2 

3.93e-06 

3.61 

2.45e-04 

3.16 

3.04e-05 

2.73 

1.23e-03 

3.71 

11 x 11 x 2 

2.76e-06 

3.71 

1.79e-04 

3.30 

2.34e-05 

2.74 

l.lle-03 

1.13 

12 x 12 x 2 

2.01e-06 

3.66 

1.33e-04 

3.35 

1.82e-05 

2.91 

8.80e-04 

2.65 

13 x 13 x 2 

1.50e-06 

3.64 

9.99e-05 

3.62 

1.43e-05 

2.97 

7.65e-04 

1.76 

14 x 14 x 2 

1.14e-06 

3.68 

7.53e-05 

3.82 

1.16e-05 

2.87 

6.45e-04 

2.30 

15 x 15 x 2 

8.81e-07 

3.77 

5.80e-05 

3.78 

9.57e-06 

2.77 

5.79e-04 

1.55 

16 x 16 x 2 

6.91e-07 

3.76 

4.53e-05 

3.81 

8.01e-06 

2.75 

4.79e-04 

2.96 


Table 10 Error convergence is shown for the fully-staggered plg = 4, plgl = 5 and conven- 
tional LGL [1,2] plgl = 4 SSDC algorithms for the isentropic Euler vortex propagation on highly 
nonuniform grids. 



Fully-staggered, 

Plg = 4, plgl 

= 5 

Conventional, plgl = 

4 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

6.59e-04 

- 

8.13e-03 

- 

2.02e-03 

- 

4.01e-02 

- 

2x2x2 

1.08e-04 

2.61 

1.98e-03 

2.04 

2.90e-04 

2.80 

5.48e-03 

2.87 

3x3x2 

5.77e-05 

1.54 

2.05e-03 

-0.09 

1.70e-04 

1.32 

8.78e-03 

-1.16 

4x4x2 

2.09e-05 

3.52 

1.40e-03 

1.32 

8.33e-05 

2.48 

4.12e-03 

2.63 

5x5x2 

5.91e-06 

5.67 

3.26e-04 

6.54 

3.27e-05 

4.19 

2.14e-03 

2.93 

6x6x2 

2.52e-06 

4.67 

7.79e-05 

7.85 

1.18e-05 

5.58 

9.91e-04 

4.23 

7x7x2 

1.23e-06 

4.63 

6.04e-05 

1.66 

8.90e-06 

1.85 

6.95e-04 

2.30 

8x8x2 

6.58e-07 

4.71 

4.17e-05 

2.77 

5.62e-06 

3.44 

3.77e-04 

4.59 

9x9x2 

3.90e-07 

4.44 

2.46e-05 

4.49 

3.59e-06 

3.82 

2.19e-04 

4.59 

10 x 10 x 2 

2.43e-07 

4.48 

1.29e-05 

6.14 

2.39e-06 

3.83 

1.64e-04 

2.75 

11 x 11 x 2 

1.56e-07 

4.67 

9.68e-06 

2.99 

1.66e-06 

3.86 

1.20e-04 

3.31 

12 x 12 x 2 

1.02e-07 

4.91 

7.30e-06 

3.24 

1.20e-06 

3.71 

7.67e-05 

5.12 

13 x 13 x 2 

6.75e-08 

5.11 

4.36e-06 

6.44 

9.03e-07 

3.55 

6.29e-05 

2.47 

14 x 14 x 2 

4.65e-08 

5.03 

2.42e-06 

7.97 

6.93e-07 

3.57 

5.40e-05 

2.06 
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Table 11 Error convergence is shown for the fully-staggered plg = 5, plgl = 6 and conven- 
tional LGL [1,2] plgl = 5 SSDC algorithms for the isentropic Euler vortex propagation on highly 
nonunifornr grids. 



Fully-staggered, 

Plg = 5, plgl 

= 6 

Conventional, plgl = 

5 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

2.93e-04 

- 

5.48e-03 

- 

6.72e-04 

- 

1.35e-02 

- 

2x2x2 

2.44e-05 

3.58 

5.74e-04 

3.26 

6.53e-05 

3.36 

2.53e-03 

2.41 

3x3x2 

1.02e-05 

2.15 

4.46e-04 

0.62 

3.72e-05 

1.39 

1.67e-03 

1.02 

4x4x2 

2.65e-06 

4.70 

8.64e-05 

5.70 

1.36e-05 

3.49 

7.67e-04 

2.70 

5x5x2 

7.36e-07 

5.75 

2.93e-05 

4.85 

5.77e-06 

3.86 

3.68e-04 

3.30 

6x6x2 

2.36e-07 

6.24 

1.24e-05 

4.73 

2.90e-06 

3.78 

1.83e-04 

3.83 

7x7x2 

1.00e-07 

5.54 

4.80e-06 

6.15 

1.06e-06 

6.51 

7.48e-05 

5.80 

8x8x2 

4.97e-08 

5.26 

2.57e-06 

4.66 

5.10e-07 

5.50 

3.09e-05 

6.62 

9x9x2 

2.55e-08 

5.67 

2.08e-06 

1.83 

2.98e-07 

4.56 

1.93e-05 

3.97 

10 x 10 x 2 

1.41e-08 

5.64 

1.25e-06 

4.83 

1.81e-07 

4.74 

1.03e-05 

5.98 

11 x 11 x 2 

8.02e-09 

5.90 

6.86e-07 

6.28 

1.16e-07 

4.66 

6.76e-06 

4.43 

12 x 12 x 2 

4.78e-09 

5.94 

3.84e-07 

6.68 

7.68e-08 

4.75 

5.05e-06 

3.35 

13 x 13 x 2 

2.99e-09 

5.86 

2.35e-07 

6.11 

5.21e-08 

4.84 

4.48e-06 

1.51 

14 x 14 x 2 

1.93e-09 

5.88 

1.76e-07 

3.88 

3.63e-08 

4.87 

3.62e-06 

2.87 


Table 12 Error convergence is shown for the fully-staggered plg = 10, plgl = H and conventional 
LGL [1,2] plgl = 10 SSDC algorithms for the isentropic Euler vortex propagation on highly 
nonunifornr grids. 



Fully-staggered, plg = 10, Plgl = 11 

Conventional, plgl = 

10 

Resolution 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

1x1x2 

1.73e-06 

- 

2.70e-05 

- 

2.87e-06 

- 

1.16e-04 

- 

2x2x2 

1.23e-08 

7.13 

4.74e-07 

5.83 

6.50e-08 

5.46 

4.50e-06 

4.69 

3x3x2 

2.71e-09 

3.73 

1.57e-07 

2.72 

1.41e-08 

3.78 

1.12e-06 

3.43 

4x4x2 

3.12e-10 

7.51 

1.49e-08 

8.19 

1.46e-09 

7.88 

1.19e-07 

7.79 

5x5x2 

3.52e-ll 

9.78 

4.26e-09 

5.62 

2.16e-10 

8.56 

2.64e-08 

6.76 

6x6x2 

4.19e-12 

11.66 

4.20e-10 

12.71 

2.99e-ll 

10.84 

2.80e-09 

12.30 
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8.2 Viscous shock 


The second test case considered in this section is the propagation of a viscous shock for which 
an exact time-dependent solution is known. The compressible Navier-Stokes equations support 
an exact solution for the viscous shock profile, under the assumption that the Prandtl number is 
Pr = |. Mass and total enthalpy are constant across a shock. Furthermore, if Pr = | then the 
momentum and energy equations are redundant. If we assume that the shock is propagating along 
the x\ coordinate direction , 16 the single momentum equation across the shock is given by 

av mk ~ (v - l)(v - v f ) = 0 ; -oo<xi<oo , t> 0; 

_ U 1 . _ U 1 , right . _ 7— 1 fJ. I'* 3 ) 

«i Mft ' f «i ,left ’ 27 Pr rh ’ 

where m is the constant mass flow across the shock. An exact solution is obtained by solving the 
momentum equation (76) for the velocity profile, v: 


x\ 


( l °9 |(v - l)(u ~Vf) | + ^jLog 



(77) 


A moving shock is recovered by applying a uniform translation to the solution. A full derivation of 
this solution appears in the thesis of Fisher [22], In this study, the values Uoo = A/ooCoo, ^oo = 2.5, 
Re oo = 10, and 7 = 1.4 are used. The viscous shock, which at t = 0 is located at the center of the 
computational domain, is propagated in the direction parallel to the horizontal axis. The domain 
is described by 

x\ € (—10, 10), X 2 € (— 10, 10), x% £ ( — 1 , 1), t > 0. 

The boundary conditions on the faces perpendicular to the horizontal axis are prescribed by penal- 
izing the numerical solution against the exact solution; periodic boundary conditions are used on 
the remaining four boundary faces of the computational domain. This is an excellent test problem 
for verifying the accuracy and functionality of the inviscid and viscous components of a compressible 
Navier-Stokes solver. 

The accuracy of the following entropy stable schemes is investigated using uniform Cartesian 
and unstructured nonuniform grids: 


• Fully-staggered SSDC algorithm with plg = 1,2, 3, 4, 5, 10 and plgl = 2, 3, 4, 5, 6 , 11. 

• Conventional LGL SSDC algorithm [1,13] with plgl = 1,2, 3, 4, 5, 10. 


8.2.1 Uniform Cartesian grids 

Different grid resolutions are examined, and the viscous shock has been propagated halfway out of 
the domain when the error measure is evaluated. This measures the effect of the penalty boundary 
condition and the interior scheme. Tables 13-18 show the convergence study for a sequence of 
maximum nine nested grids (the number of grid cells is doubled in each direction every time) . The 
number of cells in each coordinate direction is indicated in the first column of these tables (i.e. , 
“Resolution”). The L 2 norm of the error decay asymptotes towards the designed rate in each case, 
for both fully-staggered and the conventional LGL SSDC schemes. As for the Euler vortex test case, 
we note that the fully-staggered approach is more accurate than the conventional algorithm for the 

16 This is assumption is just used to simplify the notation. 
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same grid resolution. However, simple counting arguments based on inviscid flux and viscous flux 
evaluations, indicate that the cost of the staggered algorithm for a solution polynomial order p is 
comparable to that of an LGL operator [1,2], with a solution polynomial order of (p + 1). 


Table 13 Error convergence is shown for the fully-staggered plg = 1> Plgl = 2 and the conven- 
tional LGL [1,2] plgl = 1 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, 

Plg = 1 , Plgl 

= 2 

Conventional, plgl = 

1 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.61e-01 

- 

6.15e-01 

- 

6.45e-02 

- 

2.31e-01 

- 

4x4x2 

3.16e-02 

2.35 

1.24e-01 

2.31 

7.26e-02 

-0.17 

2.67e-01 

-0.21 

8x8x2 

1.27e-02 

1.32 

5.56e-02 

1.16 

3.64e-02 

1.00 

1.77e-01 

0.59 

16 x 16 x 2 

3.23e-03 

1.97 

1.87e-02 

1.57 

1.36e-02 

1.42 

7.76e-02 

1.19 

32 x 32 x 2 

8.78e-04 

1.88 

6.20e-03 

1.60 

4.23e-03 

1.69 

2.42e-02 

1.68 

64 x 64 x 2 

2.25e-04 

1.96 

1.84e-03 

1.75 

1.12e-03 

1.91 

7.60e-03 

1.67 

128 x 128 x 2 

5.62e-05 

2.00 

4.76e-04 

1.95 

2.79e-04 

2.01 

1.97e-03 

1.95 

256 x 256 x 2 

1.38e-05 

2.02 

1.18e-04 

2.01 

6.81e-05 

2.04 

4.52e-04 

2.13 

512 x 512 x 2 

3.43e-06 

2.01 

2.85e-05 

2.05 

1.67e-05 

2.03 

1.03e-04 

2.13 


Table 14 Error convergence is shown for the fully-staggered plg = 2, plgl = 3 and the conven- 
tional LGL [1,2] plgl = 2 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, 

Plg = 2, plgl 

= 3 

Conventional, plgl = 2 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

2.15e-02 

- 

9.04e-02 

- 

7.58e-02 

- 

2.08e-01 

- 

4x4x2 

8.66e-03 

1.31 

3.14e-02 

1.52 

2.29e-02 

1.73 

1.26e-01 

0.73 

8x8x2 

1.65e-03 

2.39 

8.80e-03 

1.84 

9.27e-03 

1.30 

8.48e-02 

0.57 

16 x 16 x 2 

4. Ole- 04 

2.04 

2.88e-03 

1.61 

2.62e-03 

1.82 

3.14e-02 

1.43 

32 x 32 x 2 

4.11e-05 

3.29 

2.77e-04 

3.38 

3.44e-04 

2.93 

4.15e-03 

2.92 

64 x 64 x 2 

5.33e-06 

2.95 

3.72e-05 

2.90 

4.12e-05 

3.06 

5.91e-04 

2.81 

128 x 128 x 2 

7.18e-07 

2.89 

5.21e-06 

2.83 

5.18e-06 

2.99 

7.64e-05 

2.95 

256 x 256 x 2 

1.04e-07 

2.79 

7.51e-07 

2.79 

6.69e-07 

2.95 

9.81e-06 

2.96 

512 x 512 x 2 

1.35e-08 

2.94 

1.55e-07 

2.28 

8.93e-08 

2.91 

1.28e-06 

2.94 


37 



Table 15 Error convergence is shown for the fully-staggered plg = 3, Plgl = 4 and the conven- 
tional LGL [1,2] plgl = 3 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, 

Plg = 3, plgl 

= 4 

Conventional, plgl = 

3 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

2.19e-02 

- 

8.39e-02 

- 

5.00e-02 

- 

3.15e-01 

- 

4x4x2 

3.97e-03 

2.46 

1.55e-02 

2.43 

1.13e-02 

2.15 

1.26e-01 

1.32 

8x8x2 

5.12e-04 

2.96 

3.49e-03 

2.15 

2.98e-03 

1.92 

3.43e-02 

1.87 

16 x 16 x 2 

2.65e-05 

4.27 

1.52e-04 

4.52 

1.56e-04 

4.25 

1.79e-03 

4.26 

32 x 32 x 2 

1.96e-06 

3.76 

1.38e-05 

3.46 

1.17e-05 

3.74 

1.75e-04 

3.35 

64 x 64 x 2 

1.45e-07 

3.76 

1.12e-06 

3.63 

7.79e-07 

3.91 

1.31e-05 

3.74 

128 x 128 x 2 

1.09e-08 

3.73 

8.53e-08 

3.71 

4.89e-08 

3.99 

7.85e-07 

4.07 

256 x 256 x 2 

7.08e-10 

3.95 

9.73e-09 

3.13 

3.03e-09 

4.01 

4.70e-08 

4.06 


Table 16 Error convergence is shown for the fully-staggered plg = 4, plgl = 5 and the conven- 
tional LGL [1,2] plgl = 4 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, 

Plg = 4, plgl 

= 5 

Conventional, plgl = 

4 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

6.27e-03 

- 

2.59e-02 

- 

3.75e-02 

- 

1.32e-01 

- 

4x4x2 

6.72e-04 

3.22 

4.10e-03 

2.66 

6.14e-03 

2.61 

6.88e-02 

0.94 

8x8x2 

1.33e-04 

2.34 

8.25e-04 

2.31 

3.86e-04 

3.99 

5.21e-03 

3.72 

16 x 16 x 2 

5.96e-06 

4.48 

3.87e-05 

4.41 

4.68e-05 

3.04 

7.04e-04 

2.89 

32 x 32 x 2 

1.58e-07 

5.23 

1.10e-06 

5.14 

9.90e-07 

5.56 

2.23e-05 

4.98 

64 x 64 x 2 

5.07e-09 

4.96 

4.38e-08 

4.65 

2.90e-08 

5.09 

7.67e-07 

4.86 

128 x 128 x 2 

1.69e-10 

4.91 

1.40e-09 

4.96 

9.32e-10 

4.96 

2.38e-08 

5.01 

256 x 256 x 2 

6.23e-12 

4.76 

5.69e-ll 

4.62 

3.09e-ll 

4.91 

1.07e-09 

4.48 


Table 17 Error convergence is shown for the fully-staggered plg = 5, plgl = 6 and the conven- 
tional LGL [1,2] plgl = 5 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, 

Plg = 5, plgl 

= 6 

Conventional, plgl = 

5 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

7.63e-03 

- 

3.14e-02 

- 

2.81e-02 

- 

2.16e-01 

- 

4x4x2 

6.93e-04 

3.46 

4.29e-03 

2.87 

1.69e-03 

4.06 

7.02e-03 

4.94 

8x8x2 

3.73e-05 

4.22 

2.11e-04 

4.35 

2.18e-04 

2.96 

3.29e-03 

1.10 

16 x 16 x 2 

3.30e-07 

6.82 

2.74e-06 

6.26 

1.79e-06 

6.92 

2.85e-05 

6.85 

32 x 32 x 2 

7.59e-09 

5.44 

6.09e-08 

5.49 

4.39e-08 

5.35 

1.15e-06 

4.63 

64 x 64 x 2 

1.29e-10 

5.87 

1.22e-09 

5.64 

6.70e-10 

6.03 

1.80e-08 

5.99 

128 x 128 x 2 

2.08e-12 

5.96 

2.40e-ll 

5.66 

1.08e-ll 

5.96 

3.35e-10 

5.75 
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Table 18 Error convergence is shown for the fully-staggered plg = 10, plgl = H and the conven- 
tional LGL [1,2] plgl = 10 SSDC algorithms for the viscous shock on uniform Cartesian grids. 



Fully-staggered, p LG = 10, plgl = H 

Conventional, plgl = 

10 

Resolution 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.69e-03 

- 

7.22e-03 

- 

3.05e-03 

- 

1.46e-02 

- 

4x4x2 

5.99e-06 

8.14 

3.34e-05 

7.76 

3.70e-05 

6.36 

5.75e-04 

4.67 

8x8x2 

3.95e-08 

7.25 

2.48e-07 

7.07 

1.27e-07 

8.19 

2.29e-06 

7.97 

16 x 16 x 2 

1.94e-ll 

10.99 

1.48e-10 

10.72 

6.58e-ll 

10.91 

1.64e-09 

10.45 

32 x 32 x 2 

4.96e-13 

5.29 

7.91e-12 

4.22 

4.77e-13 

7.11 

1.70e-ll 

6.58 
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8.2.2 Unstructured nonuniform grids 


A family of nested grids is constructed for the grid convergence study, by replicating N times the 
“unstructured grid kernel” shown in Figure 6. The viscous shock test case is again used to assess 
the accuracy of the fully-staggered approach on fully nonuniform grids. Tables 19-24 show the 
convergence study for a sequence of maximum six nested grids. As was the case in the previous study 
performed with the inviscid propagation of the Euler vortex, the fully-staggered SSDC approach is 
more accurate than the conventional SSDC algorithm for the same grid resolution. Furthermore, 
for the fully-staggered algorithm, the L 2 norm of the error decay asymptotes towards the designed 
rate (i.e. , (p + 1)) whereas the conventional path converges to order p. 



Figure 6. “Unstructured grid kernel” used to construct a sequence of nested grids for the viscous 
shock test case. 


Table 19 Error convergence is shown for the fully-staggered plg = 1 ) Plgl = 2 and the conven- 
tional LGL [1,2] plgl = 1 SSDC algorithms for the viscous shock on highly nonuniform unstruc- 
tured grids. 



Fully-staggered, 

Plg = 1 ) Plgl = 2 

Conventional, plgl = 

1 

Resolution 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.89e-02 

- 

1.43e-01 

- 

5.07e-02 

- 

2.56e-01 

- 

4x4x2 

7.85e-03 

1.27 

6.85e-02 

1.06 

2.57e-02 

0.98 

2.05e-01 

0.32 

8x8x2 

2.45e-03 

1.68 

3.08e-02 

1.15 

1.16e-02 

1.15 

1.04e-01 

0.98 

16 x 16 x 2 

6.46e-04 

1.92 

8.35e-03 

1.88 

6.18e-03 

0.91 

7.80e-02 

0.42 

32 x 32 x 2 

1.58e-04 

2.03 

2.15e-03 

1.96 

3.04e-03 

1.02 

4.12e-02 

0.92 

64 x 64 x 2 

4.04e-05 

1.97 

6.39e-04 

1.75 

1.92e-03 

0.67 

2.77e-02 

0.57 
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Table 20 Error convergence is shown for the fully-staggered plg = 2, plgl = 3 and the conven- 
tional LGL [1,2] plgl = 2 SSDC algorithms for the viscous shock on highly nonunifornr unstruc- 
tured grids. 



Fully-staggered, 

Plg = 2, plgl 

= 3 

Conventional, plgl = 

2 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

5.23e-03 

- 

4.10e-02 

- 

1.28e-02 

- 

1.44e-01 

- 

4x4x2 

1.03e-03 

2.34 

1.51e-02 

1.44 

4.67e-03 

1.46 

9.41e-02 

0.62 

8x8x2 

1.54e-04 

2.74 

2.67e-03 

2.49 

1.23e-03 

1.92 

3.53e-02 

1.41 

16 x 16 x 2 

2.34e-05 

2.72 

5.03e-04 

2.41 

2.02e-04 

2.61 

5.88e-03 

2.59 

32 x 32 x 2 

3.74e-06 

2.65 

8.29e-05 

2.60 

3.63e-05 

2.48 

7.34e-04 

3.00 

64 x 64 x 2 

5.30e-07 

2.82 

1.25e-05 

2.73 

8.31e-06 

2.13 

1.37e-04 

2.42 


Table 21 Error convergence is shown for the fully-staggered plg = 3, plgl = 4 and the conven- 
tional LGL [1,2] plgl = 3 SSDC algorithms for the viscous shock on highly nonuniform unstruc- 
tured grids. 



Fully-staggered, 

Plg = 3, plgl 

= 4 

Conventional, plgl = 

3 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

1.60e-03 

- 

1.74e-02 

- 

5.49e-03 

- 

1.15e-01 

- 

4x4x2 

1.79e-04 

3.16 

2.96e-03 

2.56 

1.21e-03 

2.19 

3.80e-02 

1.60 

8x8x2 

1.15e-05 

3.96 

1.42e-04 

4.39 

7.51e-05 

4.01 

2.03e-03 

4.23 

16 x 16 x 2 

8.67e-07 

3.73 

1.38e-05 

3.36 

6.23e-06 

3.59 

3.13e-04 

2.69 

32 x 32 x 2 

6.45e-08 

3.75 

1.25e-06 

3.47 

6.13e-07 

3.34 

2.90e-05 

3.44 

64 x 64 x 2 

4.65e-09 

3.79 

1.01e-07 

3.63 

7.13e-08 

3.11 

2.71e-06 

3.42 


Table 22 Error convergence is shown for the fully-staggered plg = 4, plgl = 5 and the conven- 
tional LGL [1,2] plgl = 4 SSDC algorithms for the viscous shock on highly nonuniform unstruc- 
tured grids. 



Fully-staggered, 

PLG = 4, plgl 

= 5 

Conventional, plgl = 

4 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

3.46e-04 

- 

3.82e-03 

- 

2.35e-03 

- 

6.95e-02 

- 

4x4x2 

4.68e-05 

2.88 

7.66e-04 

2.32 

1.52e-04 

3.95 

4.63e-03 

3.91 

8x8x2 

1.76e-06 

4.73 

3.83e-05 

4.32 

1.22e-05 

3.63 

6.35e-04 

2.87 

16 x 16 x 2 

6.61e-08 

4.74 

1.97e-06 

4.28 

3.89e-07 

4.98 

1.34e-05 

5.57 

32 x 32 x 2 

2.22e-09 

4.90 

7.09e-08 

4.79 

1.53e-08 

4.66 

6.96e-07 

4.26 

64 x 64 x 2 

7.45e-ll 

4.80 

3.23e-09 

4.46 

7.95e-10 

4.27 

4.03e-08 

4.11 
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Table 23 Error convergence is shown for the fully-staggered plg = 5, plgl = 6 and the conven- 
tional LGL [1,2] plgl = 5 SSDC algorithms for the viscous shock on highly nonuniform unstruc- 
tured grids. 



Fully-staggered, 

Plg = 5, plgl 

= 6 

Conventional, plgl = 5 

Resolution 

L 2 error 

L 2 rate 

L°° error L 

00 rate 

L 2 error 

L 2 rate 

L°° error 

L°° rate 

2x2x2 

2.28e-04 

- 

3.77e-03 

- 

6.79e-04 

- 

9.36e-03 

- 

4x4x2 

8.23e-06 

4.79 

2.03e-04 

4.22 

5.46e-05 

3.64 

3.18e-03 

1.56 

8x8x2 

2.25e-07 

5.19 

1.05e-05 

4.26 

5.13e-07 

6.73 

2.74e-05 

6.86 

16 x 16 x 2 

4.02e-09 

5.81 

2.32e-07 

5.51 

1.67e-08 

4.94 

1.17e-06 

4.55 

32 x 32 x 2 

6.37e-ll 

5.98 

5.05e-09 

5.52 

3.52e-10 

5.57 

3.10e-08 

5.24 


Table 24 Error convergence is shown for the fully-staggered plg = 10, Plgl = H and the con- 
ventional LGL [1, 2] plgl = 10 SSDC algorithms for the viscous shock on highly nonuniform 
unstructured grids. 



Fully-staggered, plg = 10, Plgl = H 

Conventional, plgl = 10 

Resolution 

L 2 error L 2 rate 

L°° error L°° rate 

L 2 error L 2 rate 

L°° error L°° rate 

2x2x2 

4x4x2 

8x8x2 

1.43e-06 

9.84e-09 7.18 

7.05e-12 10.45 

2.92e-05 

2.71e-07 6.75 

4.24e-10 9.32 

8.30e-06 

3.70e-08 7.81 

2.40e-ll 10.59 

5.38e-04 

2.70e-06 7.64 

2.89e-09 9.87 
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8.3 Taylor-Green vortex 


The goal of this section is to demonstrate that nonlinearly stable schemes do not require stabi- 
lization techniques (e.g., de- aliasing, filtering, limiting, over-integration procedures) for successful 
computations of under-resolved turbulent flows. 

The Taylor-Green vortex test case is used as a benchmark model problem, and is solved on the 
periodic cube [— irL < x,y,z < +ttL\. The initial condition is given by the following analytical 
expressions 

%2\ 

-jcos 



Ui = Vo sin ( — 


xi 


L 


U 2 = -To cos — sin — cos — 


X2 
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X3 
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u 3 = 0, 
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cos 


2X3 
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+ 2 


where u±, U2, u 3 are the components of the velocity in the xq-, X2-, and X 3 -directions, p is the 
pressure, and p is the density. The flow is initialized to be isothermal, i.e. , p/p = po/po = RTo- 
The Reynolds number for this flow is defined as Re = (poVqL) / p, where p is the dynamic vis- 
cosity. Starting from the initial condition, the nonlinear interactions of different flow scales yield 
vortex breakdowns. This nonlinear process is initially anisotropic and laminar; but subsequently, it 
develops into fully anisotropic turbulence that decays with the typical spectral energy distribution. 

Because our framework solves the compressible Navier-Stokes equations, we choose a Mach 
number of M = 0.08, which leads to a flow field that is essentially incompressible. This allows 
for a fair comparison with some of the incompressible simulations reported in literature. Herein, 
two Reynolds numbers are considered: Re = 800 and Re = 1,600. The Prandtl number is set to 
Pr = 0.71. 


8.3.1 Re = 800 

In this section, we report the results for Re = 800 on a uniform Cartesian grid with four hexahedrons 
in each coordinate direction. Therefore, the total number of elements in the grid is also N e = 4 3 = 
64. We compute the numerical solution with the sixteenth- (pc = 15, plgl = 16), and seventeenth- 
order ( pc = 16, plgl = 17) accurate fully-staggered algorithms. The numerical solutions presented 
herein are compared with the direct numerical simulation (DNS) of Brachet et al. [34]. Figure 
7 shows the kinetic energy dissipation rate of our computations and the reference data. The 
computation with a formally seventeenth-order accurate scheme compares very well with the DNS 
results, even on this very coarse grid. 

8.3.2 Re = 1, 600 

The simulation is run using a fully unstructured grid that contains 42 hexahedrons. The distribution 
of the element in the 2D plane is shown in Figure 8. 

Figure 9 shows the kinetic energy dissipation rate of our computations and the reference data of 
Carton de Wiart et al. [35] . The grid is too coarse to accurately resolve the flow field, however, the 
computation with a formally seventeenth-order accurate scheme is stable through all the simulation. 
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Taylor-Green vortex, Re= 800, M= 0.08 



Figure 7. Evolution of the time derivative of the kinetic energy for the Taylor-Green vortex at 
Re = 800, M = 0.08; fully-staggered SSDC algorithm. 



Figure 8. Plane distribution of the elements used for the Taylor-Green vortex at Re = 1,600, 
M = 0.08. 


This is a feat unattainable with alternative approaches based on high-order accurate linear stable 
schemes. 


8.4 Computation of a square cylinder in a supersonic free stream 

The development of a high-order accurate entropy stable discretization aims to provide the next 
generation of robust high-fidelity numerical solvers for complex fluid flow simulations, for which 
standard suboptimal algorithms suffer greatly or fail completely. By computing the flow past a 3D 
square cylinder at Re oo = 10 4 and M 0 0 = 1.5, we provide numerical evidence of such robustness for 
the complete entropy stable high-order spatial discretization. This supersonic flow is characterized 
by a very large range of length scales, strong shocks, and expansion regions that interact with each 
other, leading to complex flow patterns. During the past three decades, this fluid flow problem has 
been thoroughly investigated by several researchers for aerodynamic applications (see for instance, 
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Taylor-Green vortex, Re =1,600, M= 0.08 



Figure 9. Evolution of the time derivative of the kinetic energy for the Taylor-Green vortex at 
Re = 1,600, M = 0.08; fully-staggered SSDC algorithm. 


references 36,37). 

The domain of interest spans one square cylinder edge in the .T 3 direction, and at the two planes 
perpendicular to this coordinate direction, periodic boundary conditions are used. The flow is 
computed using an unstructured grid with 43, 936 hexahedrons. A fourth-order accurate (plg = 3, 
Plgl = 4) fully staggered entropy stable discretization without any stabilization technique is used. 
The body surface is considered adiabatic and the wall boundary conditions are imposed using the 
entropy stable approach presented in reference 2. The solution is initialized using a uniform flow 
at Moo = 1-5 with zero angle of attack. This test case was also used in reference 2 to assess the 
robustness of the conventional SSDC algorithm. 

A strong shock is formed in front of the bluff body in the beginning of the simulation. Subse- 
quently, the discontinuity moves upstream until it reaches a “stationary” position that is about 2.15 
square cylinder edges from the frontal surface of the body. During this phase, additional weaker 
shocks, which originate from the four sharp corners of the body, interact with the subsonic regions 
formed near the walls. This complicated flow pattern, yields the formation of shocklets in the wake 
of the square cylinder. 

A global view of the “high order grid,” the Mach number, density, temperature and entropy 
contours at t = 100 are shown in Figure 10. The shock has already reached a stationary position 
at t = 100 , and the flow past the square cylinder is completely unsteady, characterized by subsonic 
and supersonic regions. The formation of shocklets in the near wake region are clearly visible. 

An extensive parametric study of Mach numbers (1.1 < M < 1.8) and grids (10/F — 40A" el- 
ements) was performed using the supersonic square cylinder. All high Mach number simulations 
performed to date suggest that the staggered and conventional approaches have comparable ro- 
bustness. The failure mode for both approaches is a negative density that develops in the vicinity 
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(a) High order grid. 



(b) Mach number; AM = 0.0095. 



(c) Density; A p = 0.0090. 
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(d) Temperature; AT = 0.0077. 



(e) Entropy; As = 0.0044. 


Figure 10. Unsteady flow past a 3D square cylinder at Re oo = 10 4 and = 1.5; fourth-order 
accurate fully-staggered SSDC method (plg = 3, plgl = 4) without stabilization technique; 
t = 100. 
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of a strong shock. 


9 Conclusions 

A summation-by-parts, simultaneous-approximation-term (SBP-SAT) framework is used to de- 
velop a generalized entropy stable spectral element formulation that includes a broader selection 
of collocation points. A necessary condition for an entropy-stable, staggered operator is a restric- 
tion/prolongation interpolation pair that satisfies a precise Galerkin constraint (see equation (16)). 
This constraint is automatically satisfied when interpolating between the LG and LGL points, 
provided the polynomial order of the LGL points exceeds that of the LG points (as well as other 
less desirable combinations). Thus, the primary goals herein are to 1) extend the entropy stable 
mechanics to a staggered configuration that collocates the solution variables at the Legendre-Gauss 
(LG) points and the fluxes at the Legendre-Gauss-Lobatto (LGL) points, and 2) to establish the 
competitive advantages (if any) of the new entropy stable staggered algorithm relative to the algo- 
rithms presented in references 1,2. 

Conventional energy analysis as well as entropy analysis are used on the ID viscous Burg- 
ers’ equation, to prove the nonlinearly stability of the staggered operators for arbitrary order, 
diagonal-norm SBP operators. The comparison of the energy and entropy estimates on Burgers’ 
equation provides insight on how to proceed with the entropy stability analysis of the compressible 
Navier-Stokes equations for staggered operators. Next, the entropy analysis techniques presented 
in references 1,2 are used to develop an entropy conservative, staggered grid, spectral collocation 
operator for the 3D compressible Navier-Stokes equation. Entropy conservative/stable inviscid in- 
terface operators are then developed for the staggered formulation. The viscous interface coupling is 
based on an LDG/IP approach analogous to the coupling operators presented in references 1,2. The 
resulting spectral collocation operators are design order accurate for arbitrary order, conservative 
on the LGL points, and satisfy the additional secondary constraint of entropy stability. Extension 
to curvilinear meshes [22,23] and entropy stable solid wall BCs [13,31] follow immediately. 

Extensive numerical tests for the 3D Euler and compressible Navier-Stokes equations reveal 
that the new staggered grid entropy stable algorithms are significantly more accurate than the 
collocated Legendre-Gauss-Lobatto operators of equivalent polynomial order, that are presented in 
references 1,2. They are however, more costly to implement. The cost of a flux evaluation for the 
staggered algorithm of solution polynomial order p is comparable to that of an Legendre-Gauss- 
Lobatto points operator [1,2] of solution polynomial order (p + 1). Preliminary studies using an 
explicit temporal integrator indicate that the increased accuracy of the staggered approach nearly 
offsets the additional cost. Further investigation is required to fully establish the cost efficacy 
of over-collocated fluxes (including the spectral difference and the ffux reconstruction operators), 
relative to conventional collocated fluxes. An ongoing investigation continues that includes 1) 
the effects of implicit temporal integrators and 2) data movement on current and next generation 
hardware; computationally intensive yet extremely accurate, low memory footprint algorithms could 
be competitive in the future. 

This work provides an essential step towards an operational entropy stable framework of ar- 
bitrary order. An obvious extension of this work is the development of entropy stable, h- and/or 
refinement operators. Both scenarios require data interpolation from adjoining interfaces onto 
a common intermediate mortar, with quadrature points that do not in general coincide. A long 
term goal includes entropy stable operators for other element types including triangles, prisms and 
tetrahedra. 
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Appendix A 


Summation- by-parts 


A.l Complementary grids 

Define on the interval — 1 < x < 1, the vector of discrete points 

x = . . . ,x N -i,x N ) T , -1 < Xl,X 2 ,...,X N -l,X N <1. (Al) 

Because the approximate solution is constructed at these points, they are referred to as solution 
points. Furthermore, define a set of intermediate points prescribing bounding control volumes about 
each solution point. These (N + 1) points are referred to as flux points as they are similar in nature 
to the control volume edges employed in the finite-volume method. The distribution of the flux 
points depends on the discretization operator. The spacing between the flux points is implicitly 
contained in the SBP operator V. In fact, the diagonal elements of V are equal to the spacing 
between flux points (see figure 1), 

x = (xo,Xl, ■ . -X N ) T , X 0 =Xi, XN = XN, 

Xi Xi— i l 1? • • • i Ah 

In operator notation, this is equivalent to 

Ax = VI, 1 = (1,1,...,1) T , 

and A is as defined in equation (A6). Note that in equation (A2), the solution and flux points 
coincide at the left and right extremes of the domain. Thus 

7o = /(?i)> 7 jv = /(®v )• (A4) 

This duality is needed to define unique operators and is important in proving entropy stability. 


(A2) 

(A3) 


A. 2 Telescopic flux form 


All SBP derivative operators V can be manipulated into the telescopic flux form, 

f) f 

_(q) =p- 1 Qf + T p+ i = V~ 1 Af + T p+1 , 


(A5) 


where the N x (A + 1) matrix A is defined as 

/ —1 1 0 


A = 


0 -1 1 


0 
0 

V o 


0 

0 


0 0 \ 
0 0 


0 0 

-1 1 0 
0 -11/ 


(A6) 


The A operator calculates the undivided difference of the two adjacent fluxes. The existence of a 
telescopic form for all SBP operators is reiterated in the following lemma, presented without proof. 
(The original proof appears elsewhere [38].) 
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Lemma A.l. All differentiation matrices that satisfy the SBP convention given in equation (2) 
are telescoping operators in the V-norm and can be expressed as in equation (A5). 


This telescopic flux form admits a generalized SBP property. All SBP operators defined in 
equation (2) can be manipulated to transfer the action of the discrete derivative onto a test function 
with an equivalent order of approximation. The telescopic flux form defined in equation (A5) 
combined with the flux consistency condition (A4) results in a more generalized relation, 


4> T VV~ l A f = 0 1 (B - A)f = f(q N )f N - f(qi)<t>i ~ 


i-i 


where 


and 


A = 


o c 

— X 
1 

kj 

-l 

O C 

O C 

O C 


0 

O C 

O C 

O C 

O C 

O C 

0 

0 



0 

0 

ii 

0 

0 



0 

0 

0 

0 

0 

l 

-l 

0 


0 

0 

0 

0 

0 

0 

V o 

0 

0 

0 

l 

o / 


V 0 

0 

0 

0 

0 

i / 


A=(^) + 0(N~ 1 ) 

ox \ax J 


(A7) 


with 5x the local grid spacing. This is equivalent to the commonly used explanation of SBP in 
indicial form, 

N N—l 

Y & (fi ~ /»-l) = - f{qi)fl ~Y^i (<&+l “ ■ ( A8 ) 

i=l i= 1 


The action of the derivative is still moved onto the test function but at first order accuracy. Note 
that although this generalized property is used herein to construct entropy conservative fluxes, it 
is also instrumental for satisfying the Lax-Wendroff theorem [18] in weak form. 


A. 3 a-split fluxes 

Conservative or chain-rule fluxes constructed on the solution points x by using the SBP derivative 
operators, have a discretely equivalent representation on the flux points x [23] (see Figure 1). 
Consider the general flux of the form f(u) = V(u)W(u). Next, a-split the flux into a combination of 
conservative and chain-rule form and discretize with any diagonal norm SBP operator. The resulting 
a-split discrete operator has an equivalent telescoping flux representation of the form [23,38] 

aV~ 1 QVw + (1 — a) (VT’ _1 Qw + WP~ 1 Qv) = V^Af. (A9) 

(Again, V and W are diagonal matrices with the v and w vectors injected onto the matrix diago- 
nals.) The telescoping flux is constructed point-wise by using the expression 


N 


fi = 2 Y Y qek 

k=i + 1 1=1 


v(u e )w(u e ) + v(u k )w(u k ) , M v(ui)w(u k ) + v(u k )w(ue)~ 
a + (1 — a) 


1 < i < N - 1, 


2 v 2 
Jo = v(ui)w(ui), J N = v(u N )w(u N ), 


(A10) 


where the coefficient q^ k corresponds to the (£, k ) row and column in Q, respectively. 
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Appendix B 


Polynomial methods 


B.l Definitions 

Define on the interval — 1 < x < 1, the vectors of discrete point values, 

X = [xi, X 2 , ■ ■ • , XM- 1 , Xm] T , —1 < X\,X2,- ■ ■ ,XM-1,XM < 1, 

X = [x 1 ,X 2 ,...,X N -l,X N ] r , -1 < Xl,X 2 ,...,X N -l,X N <1. 


(Bl) 


Herein, the discrete points x and x are chosen to be the Legendre-Gauss (LG) points and the 
Legendre-Gauss-Lobatto (LGL) points, respectively. 

Next, define the interpolation operators that move data between x and x. Assume that an 
infinitely smooth function f(x) is defined on the interval — 1 < x < 1. Reading the function / and 
derivative df/dx at the discrete points x and x, we define the vectors 


f(x) = [f(x 0 )J(x i), . . . , f(xM-i),f(x M )}\ 
f ( x ) = [f{x o), /(xi), . . . , f(x jv-i), f(x at)] T , 


df 

* (x) = 

d f 

s (x) = 


df df df df ' 


(B2) 


df . df df df 


nT 


Define the Lagrange basis polynomials relative to the discrete points, x, as 

1 < j < M, 


t I \ n r M x ~ x k 

L j( x ) = Hfc=i~r 


kytj Xj Xk 

with a similar definition for the discrete points x 

r N x - x k 


l j (' x ) = n*=i 


k^j x j Xk 


1 < j < N. 


(B3) 


(B4) 


With a slight abuse of notation, define the vector of Lagrange basis polynomials relative to the 
vectors x and x as 


L(x;x) = (L 0 (x;x),L!(a;;x),...,L M _i(x;x),L M (x;x)) T , 
L(x;x) = (L 0 (x;x),L 1 (a:;x), . . . , LAr_i(a:; x), Lat(x; x)) T . 


(B5) 


This notation makes explicit reference to the set of collocation points from which the basis polyno- 
mials are derived. Finally, the vector of Lagrange basis polynomials can be evaluated at any set of 
points, thus creating a second-order tensor (matrix). For example, evaluating L(x;x) at the points 
x would yield 


L(x; x) = (L 0 (x;x),Li(x;x),...,Lm-i(x;x),Lm(x;x)) T , 


(B6) 
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where each element of the M vector being an N vector, e.g., 

Lj(x;x) = ■ ■ ,L j (x N ^i;SC),L j (x N ;Si)) T . (B7) 

Note that L(x;x) = Sjj, where Sij represents the Kronecker’s delta. 

B.2 Differentiation 

The interpolation polynomial /jv(x) (of order p = N — 1) that collocates f(x) at the points x is 
given by 

f(x) « /at(x) = L(x; x) T f(x). (B8) 

The objective is to construct collocation derivative operators in terms of the Lagrange basis poly- 
nomials on the interval. 

Theorem B.l. The derivative operator that exactly differentiates an arbitrary p-th order polyno- 
mial (p = N — 1) at the collocation points x is 

v = ( dij ) = • ( BQ ) 

Proof. Differentiating equation (B8) once yields the expression 

d Pfh = {f( x ^) T { (t). (BIO) 

ax \ ax J 

Evaluating equation (BIO) at a set of points (e.g., x), results in an expression of the form ^ = Vf, 
from which equation (12) follows immediately. □ 

Thus, the differentiation operator, V, for the collocation points x can be expressed in terms of 
the derivative of the Lagrange basis polynomials L(x;x). A Galerkin technique can also be used 
to derive an equivalent differentiation operator. 

Theorem B.2. The derivative operator that exactly differentiates an arbitrary p-th order polyno- 
mial (p = N — 1 ) at the collocation points x is 

V = (d^) = V Q. (Bll) 

Proof. First note that in addition to expression (BIO), the exact derivative, c ^- ) , of the function 
f(x) can be approximated by 

£(x)«^(x) = L te X) T £(X). (B12) 

In general, expressions (BIO) and (B12) are not equivalent. The Galerkin statement demands that 
the integral error between the two expressions be orthogonal to a set of Lagrange polynomials. 
Specifically that 

J L(®;x) ^L(x;x) T ^(x) - ^^(x;x)^ f(x)^ dx = 0 (B13) 
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which can be expressed in the equivalent form 




r| (S) - &(S) 

(B14) 

with 

II II 
<& <Gi 

( Pij ) = /^ 1 L(x;x)L(a;;x) T dx, 

(?«) = X 1 iL(x;5)(^(x;i)) T dx. 

(B15) 






Equation (Bll) would follow immediately if V is symmetric positive definite (SPD), and therefore 
invertible. The fact that V is symmetric follows immediately from definition (B15). Positive 
definiteness is more subtle. Following the definition of SPD, we pre- and post-nrultiply V by an 
arbitrary discrete vector ip, yielding 

ip T Vip = ip r L(x; x)L(x; x) T ip dx = ip(x) 2 dx, (B16) 

an expression, which is strictly greater than zero, unless ip is the null vector. Thus, the matrix V 
is SPD, therefore invertible, and (Bll) follows immediately. □ 

Remark. Once we have established that V is invertible, we can proven (Bll) directly by showing 
that V V = Q. This is accomplished by simplifying, after substituting the definitions of V, T> and 

Q. 

A proof that Q is nearly skew-symmetric is as follows. 

Theorem B.3. The matrix Q = /j L 1 L (*;*)(§ (x;x)) T dx is structurally of the form 

Q + Q T = B. (B17) 

Thus, by virtue of the structure of V and Q, the differentiation operator, V, is indeed an SBP 
operator defined by (2). 

Proof. Integrating by parts the definition of Q yields the expression 

Q = f\L(x;x)(^(x;x)) T dx = L(+l; x) L(+l; x) T - L(— 1; x) L(-l; x) T 

“ /m (^(^;x))L(x;5) T dx. 

All Lagrange polynomials based on the LGL collocation points vanish on the boundaries 
i,j < N. Thus, the boundary matrices reduce to the form 

L(+l; x) L(+l; x) T — L( — 1; x) L(— 1; x) T = S iN S jN - 8n Sji. 

Writing equation (B18) in indicial nomenclature leads to q VJ + q yi = 5n\ r djjy — <5ii Sj i, which is the 
desired result. □ 


(B18) 
for 1 < 
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B.2.1 Collocation and mass lumping 


A Legendre collocation operator may be constructed by approximating the integrals in equations 
(B15) by the LGL quadrature formula. Let 

V = (Vi,V2,---Vn-i,Vn) T (B19) 


be the nodes of the LGL quadrature formula (i.e., the zeroes of the polynomial d ^ Pn 1 ^ 1 x — [39]), 
and let uji, 1 < l < N be the quadrature weights. Define L(?^;x) as the vector of Lagrange basis 
polynomials evaluated at the quadrature point rji ; 


i.e., 


L()?z;x) = (Li(r 7 ;;x),L 2 (? 7 /;x), . . . ,L N (rji\x)) 


T 


Using these definitions, the mass and stiffness matrices V and Q are given by the expressions 

'P = Ez L (w x X L (w x )) T ^> Q = Ez L (w x )®w x ))"%- ( B2 °) 

Theorem B.4. The matrix V is diagonal for collocations points located at the LGL quadrature 
points, i.e., x = r\. Furthermore, the diagonal coefficients ofV are the integration weights lui, 1 < 
l < N used in the quadrature [21]. 

Proof. Recall that the Lagrange polynomials evaluated at the knot points satisfy the property 
Li{xj) = 5ij. Thus, the result follows immediately from the definition of the norm V = 

Ez L(?7z; x )(L(?7z; x )) T u;z. □ 


Remark. Replacing the full "P-norm in (B15) with a lower-order diagonal operator is sometimes 
referred to as “mass lumping”. 

Note that in general, V ^ V. The LGL formula is exact for polynomials of degree 2p— 1, where 
p = N — 1. However, L(x; x)(L(x; x)) T dx is of degree 2 p. Thus, the integration differs for 
the highest order term (i.e., 2p-th order term). Indeed, the two matrix norms differ by a rank one 
perturbation, i.e V = V + r y p 'D p eo(V p eo) T where eo = (1,0,--- , 0) T , V p is the highest derivative 
supported by the polynomial, and depends on polynomial order. 

The matrices Q and Q are equivalent. This follows from the fact of the two matrices are defined 
by the polynomials [^ L(x; x) {^{x; x)) dx that have a combined rank of 2 p — 1. Therefore, 
integration is exact when using the LGL integration formula. 

The uniqueness of the differentiation matrix V yields the expression 

V = V~ 1 Q = V- 1 Q. 

This statement does not contradict the fact that Indeed, the Sherman-Morrison formula [40] 

can be used to show that the difference (i.e., V — V^ 1 ) lies in the null space of the singular Q 
matrix. 


B.3 Interpolation 

We seek interpolation operators that take the discrete values f(x) of an arbitrary polynomial 
function /m(^)j from one set of points to another, (e.g., interpolating /zw(x) from x to x). Define 
the two sets of discrete points as follows: 

x = (xi,X2, ■ ■ ■ , x M -i,x M ) T , 

X = (xi,x 2 , ■ ■ ■ ,X N -l,X N ) T 


— 1 < Xl,X2, ■ ■ ■ , XM-1,XM < 1 , 
-1 < Xl,X 2 , ■ ■ ■ , X N —\, X N < 1 . 


(B21) 


57 



A Galerkin approach is used to build the interpolation operators. 

Theorem B.5. An interpolation pair that translates polynomial information between two sets of 
points x and x can be expressed as 

7 Zm-n, TlJ[_ N (B22) 


with 


i i l 

V = j L(x; x)L(x; x) T dx, V = J L(x; x)L(x; x) T dx, TZm-n = J L(x; x)L(x; x.) T dx. (B23) 

-l -l -l 


Proof. First define two distinct polynomial representations of the same function f(x), using La- 
grange polynomials and discrete function values. Assume that M and N are not equal. The 
polynomial representations are 

f(x) ~ f&[{x) = L(x; x) T f(x), f(x) « f N (x) = L(x;x) T f(x). (B24) 

Next, demand that the two polynomials approximations of f(x) be as close as possible in an 
integral sense. Two Galerkin statements for the minimization of the differences in polynomials are 

l 

J L(x;x) ^L(x; x) T f(x) — L(x; x) T f(x) j dx = 0, (B25) 

-l 

i 

J L(x; x) ^L(x; x) T f(x) — L(x;x) T f(x)j dx = 0, (B26) 

-l 

and enforce integral minimization of interpolation error between two distinct sets of Lagrange 
polynomials, (i.e. , those constructed with respect to x and x). 

Rearranging equation (B25) and (B26) yields 

f(x) = V~ x K M -n f(x) = T N -+ m f(x), f(x) = V- l n T M _ N f(Z) = T m ^>n f (x) (B27) 


with 


l i 

V = j L(x; x)L(x; Si) T dx, V = j L(x; x)L(x; x) T dx, 7 Zm-n 

-l -l 

We have already established the invertibility of V and V. 

Remark. Equation (B22) immediately implies 

VTn^m = . 


i 

j L(x; x)L(x; x) T dx. (B28) 

-l 

□ 


(B29) 
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B.3.1 Mass lumping LGL norms 

Entropy stability is achieved with a mass- lumped diagonal operator V (i.e., equation (B20)). Equa- 
tion (B29) plays a pivotal role in all staggered entropy proofs. Thus, an equivalent relation that 
holds for V = L(f?z; x)(L(rft; x)) T cuz must be derived. 

A Legendre collocation interpolation operators may be constructed by approximating the inte- 
grals in equations (B23) with the LG and LGL quadrature formula. 

Theorem B.6. The mass lumped expression 

VIlgl^lg = Zlg^lgl'P ( B3 °) 

with V = L(? 7 z; x)(L(r^; x )) T 07 is valid provided the number of LGL points exceeds by at least 

one the number of LG points (i.e., N > M). 

Proof. Define V N as the matrix T> raised to the power N. Substituting the relationship V = 
V + 7ArT ,Af eo(T ,Ar eo) into equation (B29) yields the expression 

VIlgl^lg = 1 T LG ^ LGL (V + lN V N e 0 (V N e 0 ) T ) = iJg^lglV + 1n1I g ^lglV N e o{V N e 0 ) T . 

(B31) 

If the vector condition (T> iV eo) Tlg^lgl = 0 holds, then (B30) is satisfied. 

Let £ be a matrix of basis vectors of polynomial degree M on the LG points. The action of the 
interpolation operator Zlc^lgl on £ is a rotation of each basis vector from the LG to LGL points. 
Note that the rotation does not change the polynomial order of any of the vectors. 

Now, recall that the vector V N e o is the undivided difference on the LGL. Thus, if N > M, 
then the undivided difference of a polynomial of degree < M will be the zero vector for each of the 
M basis vectors. Because the vectors are arbitrary, the vector condition ( V N eo ) Zl G ^l G l = 0 
must hold. □ 

Remark. For the work presented herein, M = p and N = p + 1. 
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